GESTIÓ DELS ESPAIS VERDS PER LA CONSERVACIÓ DE POL·LINITZADORS:

Carregar llibreries

#install.packages("vegan")
library(openxlsx)
library(readxl)
library(dplyr)
library(vegan)
library(MASS)
library(reshape2)
library(reshape)
library(ggplot2)
library(eulerr)
library(purrr)
library(stringr)
library(ggbiplot)
library(indicspecies)
library(purrr)
library(lme4)
library(Matrix)
library(lmerTest)

Carregar la base de dades original (corregida) i eliminar les dades de Mordèl·lids i de les zones “Llig”, “Mura”, i “TV-7041”

# Obtenir el directori de l'arxiu actual
directorio <- getwd()

dades_finals <- read_excel(file.path(directorio,"Carreteres_ITS+altres.xlsx"))
dades_finals <- dades_finals %>%
  filter(Familia != "Mordellidae" & !Zona %in% c("Llig", "Mura", "TV-7041"))
head(dades_finals)
## # A tibble: 6 × 25
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6        28 Campus Juliol         7  2021      4 NS         <NA>         
## # ℹ 17 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>
# Guardar la nova base de dades
write.xlsx(dades_finals, file.path(directorio, "dades_finals.xlsx"))

ANÀLISI DESCRIPTIVA

# Definir nova variable: "Mostreig", segons Zona, Mes i Any
dades_finals <- dades_finals %>%
  mutate(Mostreig = paste(Zona, Mes, Any, sep = "-"))

Riquesa

Riquesa total

riquesa_total<- dades_finals %>% 
  summarise(riquesa_total=n_distinct(ID))
riquesa_total
## # A tibble: 1 × 1
##   riquesa_total
##           <int>
## 1           137

Riquesa per cada mostreig i tractament

riquesa <- dades_finals %>%
  group_by(Mostreig, Tractament) %>%
  summarise(riquesa = n_distinct(ID))
## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.
riquesa
## # A tibble: 11 × 3
## # Groups:   Mostreig [6]
##    Mostreig             Tractament riquesa
##    <chr>                <chr>        <int>
##  1 Campus-Juliol-2021   NS              20
##  2 Campus-Juliol-2021   S               21
##  3 Campus-Jun-2021      NS              45
##  4 Campus-Jun-2021      S               34
##  5 Campus-juliol-2020   NS              39
##  6 Campus-juliol-2020   S               24
##  7 Campus-maig-2020     NS              45
##  8 Franqueses-Maig-2021 NS              20
##  9 Franqueses-Maig-2021 S               11
## 10 Sabadell-Juny-2021   NS              16
## 11 Sabadell-Juny-2021   S               12

Riquesa per cada zona i tractament

riquesa_zones <- dades_finals %>%
  group_by(Zona, Tractament) %>%
  summarise(riquesa = n_distinct(ID))
## `summarise()` has grouped output by 'Zona'. You can override using the
## `.groups` argument.
riquesa_zones
## # A tibble: 6 × 3
## # Groups:   Zona [3]
##   Zona       Tractament riquesa
##   <chr>      <chr>        <int>
## 1 Campus     NS              98
## 2 Campus     S               59
## 3 Franqueses NS              20
## 4 Franqueses S               11
## 5 Sabadell   NS              16
## 6 Sabadell   S               12

Riquesa per grups

# Abelles
riquesa_abelles<- dades_finals %>% 
  filter(Ordre == 'HymenopteraAbella') %>% 
  summarise(riquesa_abelles=n_distinct(ID))
riquesa_abelles
## # A tibble: 1 × 1
##   riquesa_abelles
##             <int>
## 1              80
# Coleòpters
riquesa_coleoptera<- dades_finals %>% 
  filter(Ordre == 'Coleoptera') %>% 
  summarise(riquesa_coleoptera=n_distinct(ID))
riquesa_coleoptera
## # A tibble: 1 × 1
##   riquesa_coleoptera
##                <int>
## 1                 33
# Sírfids
riquesa_sirfids<- dades_finals %>% 
  filter(Familia == 'Syrphidae') %>% 
  summarise(riquesa_sirfids=n_distinct(ID))
riquesa_sirfids
## # A tibble: 1 × 1
##   riquesa_sirfids
##             <int>
## 1               9

Abundància

Abundància total

abundancia_total<-dades_finals %>% 
  summarise(Numero_total=n())
abundancia_total
## # A tibble: 1 × 1
##   Numero_total
##          <int>
## 1         1111

Abundància per cada zona i tractament

abundancia_zona <- dades_finals %>%
  group_by(Zona, Tractament) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Zona'. You can override using the
## `.groups` argument.
abundancia_zona
## # A tibble: 6 × 3
## # Groups:   Zona [3]
##   Zona       Tractament Numero_total
##   <chr>      <chr>             <int>
## 1 Campus     NS                  649
## 2 Campus     S                   323
## 3 Franqueses NS                   45
## 4 Franqueses S                    14
## 5 Sabadell   NS                   29
## 6 Sabadell   S                    51

Abundància de cada gènere total

ab_genere <- dades_finals %>%
  group_by(Genere) %>%
  summarise(Numero_total = n())
head(ab_genere)
## # A tibble: 6 × 2
##   Genere     Numero_total
##   <chr>             <int>
## 1 Acmaeodera            6
## 2 Amegilla              1
## 3 Andrena               6
## 4 Anthaxia             29
## 5 Anthidium            10
## 6 Apis                 11

Abundància de cada ID per zona i tractament

abundaciazona_ID <- dades_finals %>%
  group_by(Zona, Tractament, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Zona', 'Tractament'. You can override
## using the `.groups` argument.
head(abundaciazona_ID)
## # A tibble: 6 × 4
## # Groups:   Zona, Tractament [1]
##   Zona   Tractament ID                     Numero_total
##   <chr>  <chr>      <chr>                         <int>
## 1 Campus NS         Acmaeodera cylindrica             5
## 2 Campus NS         Amegilla sp.                      1
## 3 Campus NS         Anthaxia (Anthaxia)               1
## 4 Campus NS         Anthaxia (Melanthaxia)            5
## 5 Campus NS         Anthaxia hungarica               10
## 6 Campus NS         Anthaxia sp.                      7

Abundancia de cada ID per cada mostreig i tractament

abundacia_ID <- dades_finals %>%
  group_by(Mostreig, Tractament, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.
head(abundacia_ID)
## # A tibble: 6 × 4
## # Groups:   Mostreig, Tractament [1]
##   Mostreig           Tractament ID                   Numero_total
##   <chr>              <chr>      <chr>                       <int>
## 1 Campus-Juliol-2021 NS         Amegilla sp.                    1
## 2 Campus-Juliol-2021 NS         Anthaxia sp.                    1
## 3 Campus-Juliol-2021 NS         Bombus pascuorum                1
## 4 Campus-Juliol-2021 NS         Ceratina cucurbitina            1
## 5 Campus-Juliol-2021 NS         Ceratina parvula                7
## 6 Campus-Juliol-2021 NS         Chrysididae                     1

PLOTS ABUNDÀNCIES PER GÈNERE

Campus

ab_abellescampus <- dades_finals %>%
  filter(Zona == 'Campus') %>%
  filter(Ordre == 'HymenopteraAbella') %>% 
  group_by(Genere, Tractament) %>%
  summarise(Numero_total = n(), .groups = "drop") 

p <-ggplot(ab_abellescampus, aes(y = reorder(Genere, Numero_total), x = Numero_total, fill = Tractament)) +
  geom_col(position = position_dodge(preserve = "single"),width = 0.75) +
  geom_label(aes(label = Numero_total,x= -5,y=ifelse(Tractament=="NS",as.numeric(reorder(Genere, Numero_total))-0.25,as.numeric(reorder(Genere, Numero_total))+0.2)),
             position = position_dodge(width = 0.3),
             show.legend = FALSE,
             size = 3.5,
             fill = ifelse(ab_abellescampus$Tractament=="NS","#81b29a","#e07a5f")) +
  scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +  # Ajusta los nombres de los tratamientos
  guides(fill = guide_legend(nrow = 2, byrow = TRUE)) +
  theme_bw() +
  theme(legend.position = "right",
        axis.text.y = element_text(size = 15, face = "italic"),
    axis.text.x = element_text(size = 15),  
    axis.title.x = element_text(size = 15), 
    axis.title.y = element_text(size = 15, face = "italic")
  ) +
  expand_limits(x = -6) +
  labs(x = "Abundància",
       y = "")

p

ggsave("barplot_campus.png",plot=p,width = 15, height = 11)

Sabadell

ab_abellessabadell <- dades_finals %>%
  filter(Zona == 'Sabadell') %>%
  filter(Ordre == 'HymenopteraAbella') %>% 
  group_by(Genere, Tractament) %>%
  summarise(Numero_total = n(), .groups = "drop") 

q <-ggplot(ab_abellessabadell, aes(y = reorder(Genere, Numero_total), x = Numero_total, fill = Tractament)) +
  geom_col(position = position_dodge(preserve = "single"),width = 0.75) +
  geom_label(aes(label = Numero_total,x= 0,y=ifelse(Tractament=="NS",as.numeric(reorder(Genere, Numero_total))-0.25,as.numeric(reorder(Genere, Numero_total))+0.2)),
             position = position_dodge(width = 0.3),
             show.legend = FALSE,
             size = 3.5,
             fill = ifelse(ab_abellessabadell$Tractament=="NS","#81b29a","#e07a5f")) +
  scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +  # Ajusta los nombres de los tratamientos
  guides(fill = guide_legend(nrow = 2, byrow = TRUE)) +
  theme_bw() +
  theme(legend.position = "right",
       axis.text.y = element_text(size = 15, face = "italic"),
    axis.text.x = element_text(size = 15),  
    axis.title.x = element_text(size = 15), 
    axis.title.y = element_text(size = 15, face = "italic")
  ) +
  labs(x = "Abundància",
       y = "")

q

ggsave("barplot_sabadell.png",plot=q,width = 15, height = 11)

Franqueses

ab_abellesfranqueses <- dades_finals %>%
  filter(Zona == 'Franqueses') %>%
  filter(Ordre == 'HymenopteraAbella') %>% 
  group_by(Genere, Tractament) %>%
  summarise(Numero_total = n(), .groups = "drop") 

r <-ggplot(ab_abellesfranqueses, aes(y = reorder(Genere, Numero_total), x = Numero_total, fill = Tractament)) +
  geom_col(position = position_dodge(preserve = "single"),width = 0.75) +
  geom_label(aes(label = Numero_total,x= 0,y=ifelse(Tractament=="NS",as.numeric(reorder(Genere, Numero_total))-0.25,as.numeric(reorder(Genere, Numero_total))+0.2)),
             position = position_dodge(width = 0.3),
             show.legend = FALSE,
             size = 3.5,
             fill = ifelse(ab_abellesfranqueses$Tractament=="NS","#81b29a","#e07a5f")) +
  scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +  # Ajusta los nombres de los tratamientos
  guides(fill = guide_legend(nrow = 2, byrow = TRUE)) +
   theme_bw() +
  theme(legend.position = "right",
        axis.text.y = element_text(size = 15, face = "italic"),
    axis.text.x = element_text(size = 15),  
    axis.title.x = element_text(size = 15), 
    axis.title.y = element_text(size = 15, face = "italic")
  ) +
  labs(x = "Abundància",
       y = "")

r

ggsave("barplot_franqueses.png",plot=r,width = 15, height = 11)

DIVERSITAT DE SHANNON I CORBES D’ACUMULACIÓ D’ESPÈCIES

Diversitat de Shannon per Zona i Tractament

# Preparar dades: Crear una taula de l'abundància de cada ID (columnes) segons la zona i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux2 <- cast(abundaciazona_ID, Zona + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux2 <- as.data.frame(taula_aux2)
taula_aux2[is.na(taula_aux2)] <- 0
# Agrupar Zona amb Tractament (en una nova variable Agrupacio2)
data_shannon_ZT <- taula_aux2 %>%
  mutate(Agrupacio2 = paste(Zona, Tractament, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Zona, -Tractament)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio2 i establir-la com a identificador de les files
rownames(data_shannon_ZT) <- data_shannon_ZT$Agrupacio2
data_shannon_ZT <- data_shannon_ZT %>% dplyr::select(-Agrupacio2) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig i Tractament
shannondiv2 <- diversity(data_shannon_ZT)
print(shannondiv2)
##     Campus_NS      Campus_S Franqueses_NS  Franqueses_S   Sabadell_NS 
##      3.612943      2.870658      2.627033      2.341994      2.382088 
##    Sabadell_S 
##      2.003690

Diversitat de Shannon per Mostreig i Tractament

# Preparar dades: Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux <- cast(abundacia_ID, Mostreig + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux <- as.data.frame(taula_aux)
taula_aux[is.na(taula_aux)] <- 0
# Agrupar Mostreig amb Tractament (en una nova variable Agrupacio)
data_shannon_MT <- taula_aux %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_MT) <- data_shannon_MT$Agrupacio
data_shannon_MT <- data_shannon_MT %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig i Tractament
shannondiv <- diversity(data_shannon_MT)
print(shannondiv)
##   Campus-juliol-2020_NS    Campus-juliol-2020_S   Campus-Juliol-2021_NS 
##                3.011325                2.404801                2.593731 
##    Campus-Juliol-2021_S      Campus-Jun-2021_NS       Campus-Jun-2021_S 
##                1.933200                3.161179                2.856531 
##     Campus-maig-2020_NS Franqueses-Maig-2021_NS  Franqueses-Maig-2021_S 
##                2.957046                2.627033                2.341994 
##   Sabadell-Juny-2021_NS    Sabadell-Juny-2021_S 
##                2.382088                2.003690

Corba d’acumulació d’espècies per mostreig (11 mostrejos)

corba <- specaccum(data_shannon_MT, permutations = 500)
## Warning in cor(x > 0): the standard deviation is zero
plot(corba)

Extrapolació de dades (quantes espècies més trobaríem si féssim més mostrejos)

specpool(data_shannon_MT)
##     Species     chao  chao.se jack1 jack1.se    jack2     boot  boot.se  n
## All     137 236.8148 31.40658   207 26.07681 250.1182 167.3704 12.88879 11
plot(poolaccum(data_shannon_MT))

En aquest mostreig s’han trobat 137 espècies (riquesa), però s’estima que hi ha altres espècies “unseen” (no detectades). Segons diferents estimacions, la riquesa total tenint en compte aquestes espècies no detectades seria 236 (chao), 207 (jackknife primer ordre), 250 (jackknife segon ordre), i 167 (bootstrap).

Diversitat de Shannon per Mostreig, Tractament i Trampa

# Preparar  dades:
#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa
abundacia_trampa <- dades_finals %>%
  group_by(Mostreig, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
abundacia_trampa
## # A tibble: 564 × 5
## # Groups:   Mostreig, Tractament, Trampa [101]
##    Mostreig           Tractament Trampa ID                         Numero_total
##    <chr>              <chr>       <dbl> <chr>                             <int>
##  1 Campus-Juliol-2021 NS              2 Lasioglossum glabriusculum            1
##  2 Campus-Juliol-2021 NS              4 Curculionidae                         1
##  3 Campus-Juliol-2021 NS              4 Lasioglossum malachurum               1
##  4 Campus-Juliol-2021 NS              4 Pompilidae                            1
##  5 Campus-Juliol-2021 NS              4 Pyronia cecilia                       1
##  6 Campus-Juliol-2021 NS              6 Lasioglossum glabriusculum            1
##  7 Campus-Juliol-2021 NS              6 Lasioglossum politum                  3
##  8 Campus-Juliol-2021 NS              6 Pompilidae                            1
##  9 Campus-Juliol-2021 NS              8 Chrysididae                           1
## 10 Campus-Juliol-2021 NS              8 Hylaeus sp.                           1
## # ℹ 554 more rows
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux <- cast(abundacia_trampa, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_aux <- as.data.frame(taula_aux)
taula_aux[is.na(taula_aux)] <- 0
# Agrupar Mostreig amb Tractament i Trampa
data_shannon_MTTR <- taula_aux %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament, -Trampa)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_MTTR) <- data_shannon_MTTR$Agrupacio
data_shannon_MTTR <- data_shannon_MTTR %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondiv <- diversity(data_shannon_MTTR)
print(shannondiv)
##   Campus-juliol-2020_NS_2   Campus-juliol-2020_NS_3   Campus-juliol-2020_NS_4 
##                 1.7917595                 1.8310205                 2.0449312 
##   Campus-juliol-2020_NS_6   Campus-juliol-2020_NS_8  Campus-juliol-2020_NS_10 
##                 1.7328680                 1.8310205                 1.0735428 
##  Campus-juliol-2020_NS_11  Campus-juliol-2020_NS_12  Campus-juliol-2020_NS_13 
##                 0.0000000                 1.7478681                 1.6434177 
##  Campus-juliol-2020_NS_14  Campus-juliol-2020_NS_15  Campus-juliol-2020_NS_16 
##                 1.1240357                 1.5498260                 1.3862944 
##  Campus-juliol-2020_NS_20    Campus-juliol-2020_S_1    Campus-juliol-2020_S_5 
##                 2.2718685                 0.6931472                 1.9072840 
##    Campus-juliol-2020_S_7    Campus-juliol-2020_S_9   Campus-juliol-2020_S_17 
##                 0.6924212                 0.0000000                 1.0986123 
##   Campus-juliol-2020_S_18   Campus-juliol-2020_S_19   Campus-juliol-2020_S_21 
##                 1.6094379                 1.5607104                 1.1189689 
##   Campus-juliol-2020_S_22   Campus-Juliol-2021_NS_2   Campus-Juliol-2021_NS_4 
##                 1.0751393                 0.0000000                 1.3862944 
##   Campus-Juliol-2021_NS_6   Campus-Juliol-2021_NS_8  Campus-Juliol-2021_NS_10 
##                 0.9502705                 1.4648164                 0.6931472 
##  Campus-Juliol-2021_NS_12  Campus-Juliol-2021_NS_14  Campus-Juliol-2021_NS_16 
##                 1.3296613                 1.3321790                 0.0000000 
##  Campus-Juliol-2021_NS_18  Campus-Juliol-2021_NS_20    Campus-Juliol-2021_S_1 
##                 1.3862944                 1.0986123                 1.0986123 
##    Campus-Juliol-2021_S_3    Campus-Juliol-2021_S_5    Campus-Juliol-2021_S_7 
##                 0.6365142                 2.0963526                 1.5000764 
##    Campus-Juliol-2021_S_9   Campus-Juliol-2021_S_11   Campus-Juliol-2021_S_13 
##                 1.1973401                 1.4388830                 0.0000000 
##   Campus-Juliol-2021_S_15   Campus-Juliol-2021_S_17      Campus-Jun-2021_NS_2 
##                 0.6931472                 0.5091373                 2.3345491 
##      Campus-Jun-2021_NS_4      Campus-Jun-2021_NS_6      Campus-Jun-2021_NS_8 
##                 2.0854684                 1.1537419                 2.6197294 
##     Campus-Jun-2021_NS_10     Campus-Jun-2021_NS_12     Campus-Jun-2021_NS_14 
##                 1.8310205                 2.5575077                 1.7623137 
##     Campus-Jun-2021_NS_16     Campus-Jun-2021_NS_18     Campus-Jun-2021_NS_20 
##                 1.3593685                 1.9061547                 0.9502705 
##       Campus-Jun-2021_S_1       Campus-Jun-2021_S_3       Campus-Jun-2021_S_5 
##                 2.2718685                 2.0028830                 1.7351265 
##       Campus-Jun-2021_S_7       Campus-Jun-2021_S_9      Campus-Jun-2021_S_11 
##                 0.9502705                 2.1383972                 1.5595812 
##      Campus-Jun-2021_S_13      Campus-Jun-2021_S_15      Campus-Jun-2021_S_17 
##                 1.8576598                 0.7356219                 2.1639557 
##      Campus-Jun-2021_S_19     Campus-maig-2020_NS_1     Campus-maig-2020_NS_2 
##                 1.9722470                 0.8392696                 1.3642812 
##     Campus-maig-2020_NS_3     Campus-maig-2020_NS_4     Campus-maig-2020_NS_5 
##                 1.6957425                 2.5521719                 1.5171064 
##     Campus-maig-2020_NS_6     Campus-maig-2020_NS_7     Campus-maig-2020_NS_8 
##                 2.2614951                 1.2424533                 1.7480673 
##     Campus-maig-2020_NS_9    Campus-maig-2020_NS_10    Campus-maig-2020_NS_11 
##                 1.6769878                 1.9722470                 1.7480673 
##    Campus-maig-2020_NS_12    Campus-maig-2020_NS_13    Campus-maig-2020_NS_14 
##                 1.0986123                 1.5481155                 1.8891592 
##    Campus-maig-2020_NS_15    Campus-maig-2020_NS_16    Campus-maig-2020_NS_17 
##                 1.8392967                 1.0356097                 1.3296613 
##    Campus-maig-2020_NS_18    Campus-maig-2020_NS_19    Campus-maig-2020_NS_20 
##                 1.5810938                 1.3535914                 0.9002561 
## Franqueses-Maig-2021_NS_1 Franqueses-Maig-2021_NS_2 Franqueses-Maig-2021_NS_3 
##                 1.0789922                 0.8675632                 0.6931472 
## Franqueses-Maig-2021_NS_4 Franqueses-Maig-2021_NS_5  Franqueses-Maig-2021_S_1 
##                 2.2739657                 0.7549968                 1.3321790 
##  Franqueses-Maig-2021_S_2  Franqueses-Maig-2021_S_4   Sabadell-Juny-2021_NS_1 
##                 1.0986123                 1.5607104                 1.7478681 
##   Sabadell-Juny-2021_NS_4   Sabadell-Juny-2021_NS_7   Sabadell-Juny-2021_NS_8 
##                 0.9502705                 1.7917595                 1.4750763 
##  Sabadell-Juny-2021_NS_11    Sabadell-Juny-2021_S_2    Sabadell-Juny-2021_S_3 
##                 1.3862944                 0.6365142                 0.5623351 
##    Sabadell-Juny-2021_S_5    Sabadell-Juny-2021_S_6    Sabadell-Juny-2021_S_9 
##                 1.4995094                 1.0986123                 0.6931472 
##   Sabadell-Juny-2021_S_10   Sabadell-Juny-2021_S_12 
##                 0.5004024                 1.2852930

Corba d’acumulació d’espècies segons mostreig, tractament i trampa

corba <- specaccum(data_shannon_MTTR, permutations = 500)
plot(corba)

Extrapolació de dades (quantes espècies més trobaríem si féssim més mostrejos)

specpool(data_shannon_MTTR)
##     Species     chao  chao.se    jack1 jack1.se    jack2     boot  boot.se   n
## All     137 224.1493 29.51259 201.3564 11.27322 241.7798 164.5183 5.987207 101
plot(poolaccum(data_shannon_MTTR))

Amb aquest mostreig s’han trobat 137 espècies (riquesa), però s’estima que hi ha altres espècies “unseen” (no detectades). Segons diferents estimacions, la riquesa total tenint en compte aquestes espècies no detectades seria 224 (chao), 201 (jackknife primer ordre), 241 (jackknife segon ordre), i 164 (bootstrap).

Corba i extrapolació pel Campus:

#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades pel Campus
abundacia_trampa_campus <- dades_finals %>%
  filter(Zona == 'Campus') %>%
  group_by(Mostreig, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxcampus <- cast(abundacia_trampa_campus, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxcampus <- as.data.frame(taula_auxcampus)
taula_auxcampus[is.na(taula_auxcampus)] <- 0

# Agrupar Mostreig amb Tractament i Trampa
data_shannon_campus <- taula_auxcampus %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament, -Trampa)

# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_campus) <- data_shannon_campus$Agrupacio
data_shannon_campus <- data_shannon_campus %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivcampus <- diversity(data_shannon_campus)
print(shannondivcampus)
##  Campus-juliol-2020_NS_2  Campus-juliol-2020_NS_3  Campus-juliol-2020_NS_4 
##                1.7917595                1.8310205                2.0449312 
##  Campus-juliol-2020_NS_6  Campus-juliol-2020_NS_8 Campus-juliol-2020_NS_10 
##                1.7328680                1.8310205                1.0735428 
## Campus-juliol-2020_NS_11 Campus-juliol-2020_NS_12 Campus-juliol-2020_NS_13 
##                0.0000000                1.7478681                1.6434177 
## Campus-juliol-2020_NS_14 Campus-juliol-2020_NS_15 Campus-juliol-2020_NS_16 
##                1.1240357                1.5498260                1.3862944 
## Campus-juliol-2020_NS_20   Campus-juliol-2020_S_1   Campus-juliol-2020_S_5 
##                2.2718685                0.6931472                1.9072840 
##   Campus-juliol-2020_S_7   Campus-juliol-2020_S_9  Campus-juliol-2020_S_17 
##                0.6924212                0.0000000                1.0986123 
##  Campus-juliol-2020_S_18  Campus-juliol-2020_S_19  Campus-juliol-2020_S_21 
##                1.6094379                1.5607104                1.1189689 
##  Campus-juliol-2020_S_22  Campus-Juliol-2021_NS_2  Campus-Juliol-2021_NS_4 
##                1.0751393                0.0000000                1.3862944 
##  Campus-Juliol-2021_NS_6  Campus-Juliol-2021_NS_8 Campus-Juliol-2021_NS_10 
##                0.9502705                1.4648164                0.6931472 
## Campus-Juliol-2021_NS_12 Campus-Juliol-2021_NS_14 Campus-Juliol-2021_NS_16 
##                1.3296613                1.3321790                0.0000000 
## Campus-Juliol-2021_NS_18 Campus-Juliol-2021_NS_20   Campus-Juliol-2021_S_1 
##                1.3862944                1.0986123                1.0986123 
##   Campus-Juliol-2021_S_3   Campus-Juliol-2021_S_5   Campus-Juliol-2021_S_7 
##                0.6365142                2.0963526                1.5000764 
##   Campus-Juliol-2021_S_9  Campus-Juliol-2021_S_11  Campus-Juliol-2021_S_13 
##                1.1973401                1.4388830                0.0000000 
##  Campus-Juliol-2021_S_15  Campus-Juliol-2021_S_17     Campus-Jun-2021_NS_2 
##                0.6931472                0.5091373                2.3345491 
##     Campus-Jun-2021_NS_4     Campus-Jun-2021_NS_6     Campus-Jun-2021_NS_8 
##                2.0854684                1.1537419                2.6197294 
##    Campus-Jun-2021_NS_10    Campus-Jun-2021_NS_12    Campus-Jun-2021_NS_14 
##                1.8310205                2.5575077                1.7623137 
##    Campus-Jun-2021_NS_16    Campus-Jun-2021_NS_18    Campus-Jun-2021_NS_20 
##                1.3593685                1.9061547                0.9502705 
##      Campus-Jun-2021_S_1      Campus-Jun-2021_S_3      Campus-Jun-2021_S_5 
##                2.2718685                2.0028830                1.7351265 
##      Campus-Jun-2021_S_7      Campus-Jun-2021_S_9     Campus-Jun-2021_S_11 
##                0.9502705                2.1383972                1.5595812 
##     Campus-Jun-2021_S_13     Campus-Jun-2021_S_15     Campus-Jun-2021_S_17 
##                1.8576598                0.7356219                2.1639557 
##     Campus-Jun-2021_S_19    Campus-maig-2020_NS_1    Campus-maig-2020_NS_2 
##                1.9722470                0.8392696                1.3642812 
##    Campus-maig-2020_NS_3    Campus-maig-2020_NS_4    Campus-maig-2020_NS_5 
##                1.6957425                2.5521719                1.5171064 
##    Campus-maig-2020_NS_6    Campus-maig-2020_NS_7    Campus-maig-2020_NS_8 
##                2.2614951                1.2424533                1.7480673 
##    Campus-maig-2020_NS_9   Campus-maig-2020_NS_10   Campus-maig-2020_NS_11 
##                1.6769878                1.9722470                1.7480673 
##   Campus-maig-2020_NS_12   Campus-maig-2020_NS_13   Campus-maig-2020_NS_14 
##                1.0986123                1.5481155                1.8891592 
##   Campus-maig-2020_NS_15   Campus-maig-2020_NS_16   Campus-maig-2020_NS_17 
##                1.8392967                1.0356097                1.3296613 
##   Campus-maig-2020_NS_18   Campus-maig-2020_NS_19   Campus-maig-2020_NS_20 
##                1.5810938                1.3535914                0.9002561
corbacampus <- specaccum(data_shannon_campus, permutations = 500)
plot(corbacampus)

# Extrapolació de les dades
specpool(data_shannon_campus)
##     Species     chao  chao.se    jack1 jack1.se    jack2     boot  boot.se  n
## All     118 184.0553 24.55884 170.3457 9.341038 201.8116 140.7347 5.098803 81

Corba i extrapolació per Franqueses:

#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades per Franqueses
abundacia_trampa_franqueses <- dades_finals %>%
  filter(Zona == 'Franqueses') %>%
  group_by(Mostreig, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxfranq <- cast(abundacia_trampa_franqueses, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxfranq <- as.data.frame(taula_auxfranq)
taula_auxfranq[is.na(taula_auxfranq)] <- 0

# Agrupar Mostreig amb Tractament i Trampa
data_shannon_franqueses <- taula_auxfranq %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament, -Trampa)

# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_franqueses) <- data_shannon_franqueses$Agrupacio
data_shannon_franqueses <- data_shannon_franqueses %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivfranq <- diversity(data_shannon_franqueses)
print(shannondivfranq)
## Franqueses-Maig-2021_NS_1 Franqueses-Maig-2021_NS_2 Franqueses-Maig-2021_NS_3 
##                 1.0789922                 0.8675632                 0.6931472 
## Franqueses-Maig-2021_NS_4 Franqueses-Maig-2021_NS_5  Franqueses-Maig-2021_S_1 
##                 2.2739657                 0.7549968                 1.3321790 
##  Franqueses-Maig-2021_S_2  Franqueses-Maig-2021_S_4 
##                 1.0986123                 1.5607104
corbafranqueses <- specaccum(data_shannon_franqueses, permutations = 500)
plot(corbafranqueses)

# Extrapolació de les dades
specpool(data_shannon_franqueses)
##     Species     chao  chao.se  jack1 jack1.se    jack2     boot  boot.se n
## All      26 52.32292 16.96108 42.625 8.605049 53.01786 33.13315 4.330973 8

Corba i extrapolació per Sabadell

#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades per Sabadell
abundacia_trampa_sabadell <- dades_finals %>%
  filter(Zona == 'Sabadell') %>%
  group_by(Mostreig, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxsaba <- cast(abundacia_trampa_sabadell, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxsaba <- as.data.frame(taula_auxsaba)
taula_auxsaba[is.na(taula_auxsaba)] <- 0

# Agrupar Mostreig amb Tractament i Trampa
data_shannon_sabadell <- taula_auxsaba %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament, -Trampa)

# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_sabadell) <- data_shannon_sabadell$Agrupacio
data_shannon_sabadell <- data_shannon_sabadell %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivsaba <- diversity(data_shannon_sabadell)
print(shannondivsaba)
##  Sabadell-Juny-2021_NS_1  Sabadell-Juny-2021_NS_4  Sabadell-Juny-2021_NS_7 
##                1.7478681                0.9502705                1.7917595 
##  Sabadell-Juny-2021_NS_8 Sabadell-Juny-2021_NS_11   Sabadell-Juny-2021_S_2 
##                1.4750763                1.3862944                0.6365142 
##   Sabadell-Juny-2021_S_3   Sabadell-Juny-2021_S_5   Sabadell-Juny-2021_S_6 
##                0.5623351                1.4995094                1.0986123 
##   Sabadell-Juny-2021_S_9  Sabadell-Juny-2021_S_10  Sabadell-Juny-2021_S_12 
##                0.6931472                0.5004024                1.2852930
corbasabadell <- specaccum(data_shannon_sabadell, permutations = 500)
plot(corbasabadell)

# Extrapolació de les dades
specpool(data_shannon_sabadell)
##     Species     chao  chao.se    jack1 jack1.se   jack2     boot  boot.se  n
## All      23 45.45833 17.10731 35.83333 5.316379 44.4697 28.44918 2.631169 12

T-student per comparar la diversitat de Shannon entre mostrejos i tractaments

#Convertir la taula de shannondiv en dataframe
shannon_ttest <- as.data.frame(shannondiv)

# Convertir l'ID de les files (mostreig + tractament + número de trampa) en una variable (nova columna)
shannon_ttest$row <- rownames(shannon_ttest)

# De la nova variable, separar en dues columnes el mostreig i eltractament (i no utilitzem el número de trampa)
shannon_ttest <- shannon_ttest %>% mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2], Mostreig = str_split(row, "_", simplify = TRUE)[,1]) %>% dplyr::select(-row)

# Prova t-student i boxplots
shannonT <- shannon_ttest %>% 
  group_by(Mostreig) %>% 
  summarise(
    mean_S = mean(shannondiv[Tractament=='S']),
    mean_NS = mean(shannondiv[Tractament=='NS']),
    t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(shannondiv ~ Tractament)$p.value,
                             NA),.groups = "drop")
  
library(openxlsx)
write.xlsx(shannonT,"shannonT.xlsx")
  
ggplot(shannon_ttest, aes(x = factor(Mostreig), y = shannondiv, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Diversitat de Shannon x Zona i Tractaments",
     x = "Mostreig",
     y = "Diversitat Shannon") +
  scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
  scale_color_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
  theme_minimal()

Només en el mostreig de Sabadell la diversitat de Shannon és significativament diferent (p valor < 0.05): l’índex de Shannon és major al tractament NS.

Diversitat de Shannon d’abelles (per mostreig)

# Preparar les dades: crear taula de dades d'abundància només per abelles
abundanciabelles_ID <- dades_finals %>% 
  mutate(Mostreig = paste(Zona, Mes, Any, sep = "-")) %>% 
  filter(Ordre == 'HymenopteraAbella') %>% 
  group_by(Mostreig, Tractament, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.
# View(abundanciabelles_ID)

#Crear una taula auxiliar de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux <- cast(abundanciabelles_ID, Mostreig + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux <- as.data.frame(taula_aux)
taula_aux[is.na(taula_aux)] <- 0

#Agrupar Mostreig i tractament en una nova variable (Agrupacio)
data_shannon_abelles <- taula_aux %>% 
  mutate(Agrupacio = paste(Mostreig, Tractament, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament)

#Definir els rownames amb la variable Agrupacio (mostreig i tractament)
rownames(data_shannon_abelles) <- data_shannon_abelles$Agrupacio
data_shannon_abelles <- data_shannon_abelles %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

#View(data_shannon_abelles)
# Càlcul diversitat de Shannon
shannondiv_abelles <- diversity(data_shannon_abelles)
print(shannondiv_abelles)
##   Campus-juliol-2020_NS    Campus-juliol-2020_S   Campus-Juliol-2021_NS 
##                2.611229                1.905899                2.090031 
##    Campus-Juliol-2021_S      Campus-Jun-2021_NS       Campus-Jun-2021_S 
##                1.806531                2.701045                2.080060 
##     Campus-maig-2020_NS Franqueses-Maig-2021_NS  Franqueses-Maig-2021_S 
##                2.199104                2.239746                1.560710 
##   Sabadell-Juny-2021_NS    Sabadell-Juny-2021_S 
##                2.045357                1.884985

Shannon abelles per trampa

#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa
abelles_abundacia_trampa <- dades_finals %>%
  filter(Ordre == 'HymenopteraAbella') %>% 
  group_by(Mostreig, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
abelles_abundacia_trampa
## # A tibble: 371 × 5
## # Groups:   Mostreig, Tractament, Trampa [100]
##    Mostreig           Tractament Trampa ID                         Numero_total
##    <chr>              <chr>       <dbl> <chr>                             <int>
##  1 Campus-Juliol-2021 NS              2 Lasioglossum glabriusculum            1
##  2 Campus-Juliol-2021 NS              4 Lasioglossum malachurum               1
##  3 Campus-Juliol-2021 NS              6 Lasioglossum glabriusculum            1
##  4 Campus-Juliol-2021 NS              6 Lasioglossum politum                  3
##  5 Campus-Juliol-2021 NS              8 Lasioglossum glabriusculum            3
##  6 Campus-Juliol-2021 NS              8 Megachile argentata                   3
##  7 Campus-Juliol-2021 NS              8 Seladonia gemmea                      1
##  8 Campus-Juliol-2021 NS             10 Bombus pascuorum                      1
##  9 Campus-Juliol-2021 NS             12 Halictus sp.                          1
## 10 Campus-Juliol-2021 NS             12 Lasioglossum glabriusculum            2
## # ℹ 361 more rows
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux_abelles <- cast(abelles_abundacia_trampa, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_aux_abelles <- as.data.frame(taula_aux_abelles)
taula_aux_abelles[is.na(taula_aux_abelles)] <- 0
# Agrupar Mostreig amb Tractament i Trampa
abelles_shannon_MTTR <- taula_aux_abelles %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament, -Trampa)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(abelles_shannon_MTTR) <- abelles_shannon_MTTR$Agrupacio
abelles_shannon_MTTR <- abelles_shannon_MTTR %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
abelles_shannondiv <- diversity(abelles_shannon_MTTR)
print(abelles_shannondiv)
##   Campus-juliol-2020_NS_2   Campus-juliol-2020_NS_3   Campus-juliol-2020_NS_4 
##                 1.3862944                 1.6674619                 1.6417347 
##   Campus-juliol-2020_NS_6   Campus-juliol-2020_NS_8  Campus-juliol-2020_NS_10 
##                 1.3296613                 1.6674619                 0.7963116 
##  Campus-juliol-2020_NS_11  Campus-juliol-2020_NS_12  Campus-juliol-2020_NS_13 
##                 0.0000000                 1.5607104                 1.6434177 
##  Campus-juliol-2020_NS_14  Campus-juliol-2020_NS_15  Campus-juliol-2020_NS_16 
##                 0.7549968                 1.5498260                 0.6931472 
##  Campus-juliol-2020_NS_20    Campus-juliol-2020_S_1    Campus-juliol-2020_S_5 
##                 1.7478681                 0.6931472                 1.3517840 
##    Campus-juliol-2020_S_7    Campus-juliol-2020_S_9   Campus-juliol-2020_S_17 
##                 0.5505734                 0.0000000                 1.0986123 
##   Campus-juliol-2020_S_18   Campus-juliol-2020_S_19   Campus-juliol-2020_S_21 
##                 1.3862944                 0.6931472                 1.1189689 
##   Campus-juliol-2020_S_22   Campus-Juliol-2021_NS_2   Campus-Juliol-2021_NS_4 
##                 0.6108643                 0.0000000                 0.0000000 
##   Campus-Juliol-2021_NS_6   Campus-Juliol-2021_NS_8  Campus-Juliol-2021_NS_10 
##                 0.5623351                 1.0042425                 0.0000000 
##  Campus-Juliol-2021_NS_12  Campus-Juliol-2021_NS_14  Campus-Juliol-2021_NS_16 
##                 1.0549202                 1.0397208                 0.0000000 
##  Campus-Juliol-2021_NS_18  Campus-Juliol-2021_NS_20    Campus-Juliol-2021_S_1 
##                 1.3862944                 0.6931472                 0.6931472 
##    Campus-Juliol-2021_S_3    Campus-Juliol-2021_S_5    Campus-Juliol-2021_S_7 
##                 0.6365142                 1.8662160                 1.4369665 
##    Campus-Juliol-2021_S_9   Campus-Juliol-2021_S_11   Campus-Juliol-2021_S_13 
##                 1.1973401                 1.4388830                 0.0000000 
##   Campus-Juliol-2021_S_15   Campus-Juliol-2021_S_17      Campus-Jun-2021_NS_2 
##                 0.6931472                 0.5091373                 1.4735024 
##      Campus-Jun-2021_NS_4      Campus-Jun-2021_NS_6      Campus-Jun-2021_NS_8 
##                 1.0986123                 0.8675632                 2.3393717 
##     Campus-Jun-2021_NS_10     Campus-Jun-2021_NS_12     Campus-Jun-2021_NS_14 
##                 1.3862944                 2.0947290                 1.4925477 
##     Campus-Jun-2021_NS_16     Campus-Jun-2021_NS_18     Campus-Jun-2021_NS_20 
##                 1.2130076                 1.5607104                 0.0000000 
##       Campus-Jun-2021_S_1       Campus-Jun-2021_S_3       Campus-Jun-2021_S_5 
##                 1.6094379                 1.5498260                 1.0397208 
##       Campus-Jun-2021_S_7       Campus-Jun-2021_S_9      Campus-Jun-2021_S_11 
##                 0.5623351                 1.5810938                 1.5595812 
##      Campus-Jun-2021_S_13      Campus-Jun-2021_S_15      Campus-Jun-2021_S_17 
##                 1.4388830                 0.4101163                 1.3321790 
##      Campus-Jun-2021_S_19     Campus-maig-2020_NS_1     Campus-maig-2020_NS_2 
##                 1.2424533                 0.6682485                 0.9611277 
##     Campus-maig-2020_NS_3     Campus-maig-2020_NS_4     Campus-maig-2020_NS_5 
##                 0.6365142                 1.6094379                 0.6365142 
##     Campus-maig-2020_NS_6     Campus-maig-2020_NS_7     Campus-maig-2020_NS_8 
##                 1.8845210                 0.6931472                 1.3862944 
##     Campus-maig-2020_NS_9    Campus-maig-2020_NS_10    Campus-maig-2020_NS_11 
##                 1.3321790                 0.9502705                 1.3862944 
##    Campus-maig-2020_NS_12    Campus-maig-2020_NS_13    Campus-maig-2020_NS_14 
##                 0.6931472                 1.3114313                 1.3296613 
##    Campus-maig-2020_NS_15    Campus-maig-2020_NS_16    Campus-maig-2020_NS_17 
##                 1.4978661                 0.9502705                 0.6365142 
##    Campus-maig-2020_NS_18    Campus-maig-2020_NS_19 Franqueses-Maig-2021_NS_1 
##                 1.3862944                 0.9251291                 0.6931472 
## Franqueses-Maig-2021_NS_2 Franqueses-Maig-2021_NS_3 Franqueses-Maig-2021_NS_4 
##                 0.5004024                 0.6931472                 1.6674619 
## Franqueses-Maig-2021_NS_5  Franqueses-Maig-2021_S_1  Franqueses-Maig-2021_S_2 
##                 1.0986123                 0.0000000                 1.0986123 
##  Franqueses-Maig-2021_S_4   Sabadell-Juny-2021_NS_1   Sabadell-Juny-2021_NS_4 
##                 0.6931472                 1.7478681                 0.0000000 
##   Sabadell-Juny-2021_NS_7   Sabadell-Juny-2021_NS_8  Sabadell-Juny-2021_NS_11 
##                 1.3862944                 1.2424533                 1.3862944 
##    Sabadell-Juny-2021_S_2    Sabadell-Juny-2021_S_3    Sabadell-Juny-2021_S_5 
##                 0.6365142                 0.5623351                 1.4995094 
##    Sabadell-Juny-2021_S_6    Sabadell-Juny-2021_S_9   Sabadell-Juny-2021_S_10 
##                 0.6931472                 0.6931472                 0.5004024 
##   Sabadell-Juny-2021_S_12 
##                 1.0986123
#Convertir la taula de shannondiv en dataframe
abelles_shannon_ttest <- as.data.frame(abelles_shannondiv)

# Convertir l'ID de les files (mostreig + tractament + número de trampa) en una variable (nova columna)
abelles_shannon_ttest$row <- rownames(abelles_shannon_ttest)

# De la nova variable, separar en dues columnes el mostreig i eltractament (i no utilitzem el número de trampa)
abelles_shannon_ttest <- abelles_shannon_ttest %>% mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2], Mostreig = str_split(row, "_", simplify = TRUE)[,1]) %>% dplyr::select(-row)

# Prova t-student i boxplots
abelles_shannonT <- abelles_shannon_ttest %>% 
  group_by(Mostreig) %>% 
  summarise(
    mean_S = mean(abelles_shannondiv[Tractament=='S']),
    mean_NS = mean(abelles_shannondiv[Tractament=='NS']),
    t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(abelles_shannondiv ~ Tractament)$p.value,
                             NA),.groups = "drop")
  
write.xlsx(abelles_shannonT,"abelles_shannonT.xlsx")
  
ggplot(abelles_shannon_ttest, aes(x = factor(Mostreig), y = abelles_shannondiv, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Diversitat de Shannon d'abelles x Zona i Tractaments",
     x = "Mostreig",
     y = "Diversitat Shannon") +
  scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
  scale_color_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
  theme_minimal()

Corbes d’acumulació d’espècies per mostreig

Corba d’acumulació. Riquesa en funció del nº de trampes (specaccum)

# A la taula data_shannon_MTTR extreure una nova columna que indiqui el mostreig
 metadata_shannon <- data_shannon_MTTR %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) 

# View(metadata_shannon)

# Subset each mostreig into its own dataframe (i eliminar la columna de "mostreig" perquè totes les dades de la taula siguin numèriques)
metadata_shannon %>% filter(Zona == 'Campus') %>% dplyr::select(-Zona)-> Campus
metadata_shannon %>% filter(Zona == 'Franqueses') %>% dplyr::select(-Zona) -> Franqueses
metadata_shannon %>% filter(Zona == 'Llig') %>% dplyr::select(-Zona) -> Llig
metadata_shannon %>% filter(Zona == 'Mura') %>% dplyr::select(-Zona) -> Mura
metadata_shannon %>% filter(Zona == 'Sabadell') %>% dplyr::select(-Zona) -> Sabadell
metadata_shannon %>% filter(Zona == 'TV') %>% dplyr::select(-Zona) -> TV7041

# Corba d'acumulació d'espècies per cada mostreig
corba1 = specaccum(Campus, method = "random")
corba2 = specaccum(Franqueses, method = "random")
corba3 = specaccum(Sabadell, method = "random")

# Creating data frames for each set of data
data1 <- data.frame(Sites = corba1$sites, Richness = corba1$richness, SD = corba1$sd)
data1$Zona <- "Campus"

data2 <- data.frame(Sites = corba2$sites, Richness = corba2$richness, SD = corba2$sd)
data2$Zona <- "Franqueses"

data3 <- data.frame(Sites = corba3$sites, Richness = corba3$richness, SD = corba3$sd)
data3$Zona <- "Sabadell"


# Combine all data frames into one
dades_combinades <- rbind(data1, data2, data3)

library(ggplot2)

# Gráfico de dispersión con líneas y error bars
ggplot(dades_combinades, aes(x = Sites, y = Richness, color = Zona)) +
  geom_point() +
  geom_line(aes(group = Zona), alpha = 0.7) +  # Líneas
  # geom_ribbon(aes(ymin = Richness - 2 * SD, ymax = Richness + 2 * SD, fill = Mostreig), alpha = 0.01) +  # Incertidumbre
  geom_errorbar(aes(ymin = Richness - 2 * SD, ymax = Richness + 2 * SD), width = 0.2, alpha = 0.5) +  # Barras de error
  ylim(0,NA)+
  labs(x = "Nº trampes", y = "Riquesa", color = "Zona") +
  theme_minimal()

Corba de rarefracció. Riquesa en funció del nombre d’individus (rarefy)

dades_agrupades <- metadata_shannon %>%
  group_by(Zona) %>%
  summarise(across(everything(), sum)) %>% 
  as.data.frame()

#View (dades_agrupades) 

# Buscar el nombre mínim d'observacions en una zona (és 84, a TV juliol)
abundacia_trampa <- dades_finals %>%
  group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.
abundacia_trampa
## # A tibble: 564 × 7
## # Groups:   Zona, Mes, Any, Tractament, Trampa [101]
##    Zona   Mes      Any Tractament Trampa ID                         Numero_total
##    <chr>  <chr>  <dbl> <chr>       <dbl> <chr>                             <int>
##  1 Campus Juliol  2021 NS              2 Lasioglossum glabriusculum            1
##  2 Campus Juliol  2021 NS              4 Curculionidae                         1
##  3 Campus Juliol  2021 NS              4 Lasioglossum malachurum               1
##  4 Campus Juliol  2021 NS              4 Pompilidae                            1
##  5 Campus Juliol  2021 NS              4 Pyronia cecilia                       1
##  6 Campus Juliol  2021 NS              6 Lasioglossum glabriusculum            1
##  7 Campus Juliol  2021 NS              6 Lasioglossum politum                  3
##  8 Campus Juliol  2021 NS              6 Pompilidae                            1
##  9 Campus Juliol  2021 NS              8 Chrysididae                           1
## 10 Campus Juliol  2021 NS              8 Hylaeus sp.                           1
## # ℹ 554 more rows
min_n <- abundacia_trampa %>% 
  group_by(Zona) %>% 
  summarize (n_especimens = sum(Numero_total)) %>% 
  summarize(min=min(n_especimens)) %>% 
  pull(min)

rownames(dades_agrupades) <- dades_agrupades$Zona
dades_agrupades <- dades_agrupades[,-1]
dades_agrupades <- dades_agrupades %>% 
                    mutate_all(as.numeric)

vegan::rarefy(dades_agrupades,sample=min_n)
##     Campus Franqueses   Sabadell 
##   28.12144   26.00000   19.50913 
## attr(,"Subsample")
## [1] 59
col <- c("#81b29a","#3d405b", "#f2cc8f")

rarecurve(dades_agrupades, step = 5, col = col, lwd = 7, cex = 0.01, cex.axis = 2, cex.lab = 2)
legend("right", legend = rownames(dades_agrupades), col = col, lwd = 5, cex = 2, bty = "n")

Corba d’interpolació: riquesa en funció del nombre d’individus (iNEXT). (??)

library(iNEXT)
transposades <- t(dades_agrupades)
D_abund <- iNEXT (transposades, datatype = 'abundance')
plot (D_abund)

Corbes de rarefracció d’abelles per cada tractament

# Extreure la columna de tractament a la taula d'abundàncies per ID de cada mostreig
dades_shannon_abelles <- data_shannon_abelles %>% 
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2])

# Agrupar per tractament
dades_agrupades_tractament <- dades_shannon_abelles %>%
  group_by(Tractament) %>%
  summarise(across(everything(), sum)) %>% 
  as.data.frame()

#View (dades_agrupades_tractament) 

dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]

# Buscar el nombre mínim d'observacions en una zona
abundacia_trampa_abelles <- dades_abelles %>%
  group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.
#View (abundacia_trampa_abelles)

min_n <- abundacia_trampa_abelles %>% 
  group_by(Tractament) %>% 
  summarize (n_especimens = sum(Numero_total)) %>% 
  summarize(min=min(n_especimens)) %>% 
  pull(min)

rownames(dades_agrupades_tractament) <- dades_agrupades_tractament$Tractament
dades_agrupades_tractament <- dades_agrupades_tractament[,-1]
dades_agrupades_tractament <- dades_agrupades_tractament %>% 
                    mutate_all(as.numeric)

vegan::rarefy(dades_agrupades_tractament,sample=min_n)
##       NS        S 
## 59.89088 37.00000 
## attr(,"Subsample")
## [1] 321
col <- c("#81b29a","#e07a5f")

rarecurve(dades_agrupades_tractament, step = 5, col = col, lwd = 2, cex = 0.01, cex.axis = 1.1, cex.lab = 1.1)
legend("right", legend = rownames(dades_agrupades_tractament), col = col, lwd = 3, cex = 1.1, bty = "n")

Corbes de rarefracció d’abelles al Campus per tractament

# Extreure la columna Zona de dades_shannon_abelles i filtrar per Campus
dades_shannon_abellescampus <- dades_shannon_abelles %>% 
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1])
dades_shannon_abellescampus <- dades_shannon_abellescampus[dades_shannon_abellescampus$Zona == "Campus", ]

# Agrupar per tractament
dadescampus_agrupades_tractament <- dades_shannon_abellescampus %>%
  group_by(Tractament) %>%
  summarise(across(where(is.numeric), sum)) %>% 
  as.data.frame()

#View (dadescampus_agrupades_tractament) 

# Buscar el nombre mínim d'observacions en una zona
abundacia_trampa_abellescampus <- dades_abelles %>%
  group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.
#View (abundacia_trampa_abellescampus)

min_n <- abundacia_trampa_abellescampus %>% 
  group_by(Tractament) %>% 
  summarize (n_especimens = sum(Numero_total)) %>% 
  summarize(min=min(n_especimens)) %>% 
  pull(min)

rownames(dadescampus_agrupades_tractament) <- dadescampus_agrupades_tractament$Tractament
dadescampus_agrupades_tractament <- dadescampus_agrupades_tractament[,-1]
dadescampus_agrupades_tractament <- dadescampus_agrupades_tractament %>% 
                    mutate_all(as.numeric)

vegan::rarefy(dadescampus_agrupades_tractament,sample=min_n)
## Warning in vegan::rarefy(dadescampus_agrupades_tractament, sample = min_n):
## requested 'sample' was larger than smallest site maximum (266)
##       NS        S 
## 54.01949 35.00000 
## attr(,"Subsample")
## [1] 321
col <- c("#81b29a","#e07a5f")
rarecurve(dadescampus_agrupades_tractament, step = 5, col = col, lwd = 4, cex = 0.01, cex.axis = 1.1, cex.lab = 1.1)
legend("right", legend = rownames(dadescampus_agrupades_tractament), col = col, lwd = 3, cex = 1.1, bty = "n")

PROVES T-STUDENT

Comparació de la riquesa per mostrejos i tractaments

Riquesa total

riquesa_x_trampa <- dades_finals %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop") 
  
  proves_t <- riquesa_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Riquesa[Tractament == "S"]),
      mean_NS = mean(Riquesa[Tractament == "NS"]),
      t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(Riquesa ~ Tractament)$p.value,
                             NA)
      ,.groups = "drop")
  
  proves_t
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     4.33    3.1         0.303 
## 2 Campus-Jun-2021        7.3     9.5         0.190 
## 3 Campus-juliol-2020     3.89    5.92        0.0281
## 4 Campus-maig-2020     NaN       6.55       NA     
## 5 Franqueses-Maig-2021   4       4.6         0.743 
## 6 Sabadell-Juny-2021     3.43    4.8         0.201
ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA, aes(color = Tractament)) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Riquesa per Zona i Tractaments",
       x = "Mostreig",
       y = "Riquesa") +
  scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
  scale_color_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
  theme_minimal()

Generalment es compleix un patró de major riquesa al tractament NS (excepte a Campus-Juliol-2021). Al mostreig del Campus-Juliol_2020 aquesta diferència és significativa (p-valor= 0,0281)

Riquesa d’abelles

# Filtrar les dades d'abelles
  dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
head(dades_abelles)  
## # A tibble: 6 × 26
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
riquesa_x_trampa <- dades_abelles %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop") 
  
  proves_t_abelles <- riquesa_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Riquesa[Tractament == "S"]),
      mean_NS = mean(Riquesa[Tractament == "NS"]),
      t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(Riquesa ~ Tractament)$p.value,
                             NA)
      ,.groups = "drop")
  
  proves_t_abelles
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     3.89    2.1        0.0847 
## 2 Campus-Jun-2021        4.3     5.4        0.338  
## 3 Campus-juliol-2020     2.78    4.46       0.00961
## 4 Campus-maig-2020     NaN       3.84      NA      
## 5 Franqueses-Maig-2021   2       3          0.341  
## 6 Sabadell-Juny-2021     3.14    3.8        0.568
ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Riquesa d'abelles per Zona i Tractaments",
       x = "Mostreig",
       y = "Riquesa") +
  scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
  scale_color_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
  theme_minimal()

Generalment es compleix un patró de major riquesa d’abelles al tractament NS (excepte a Campus-Juliol-2021). Al mostreig del Campus-Juliol_2020 aquesta diferència és significativa (p-valor= 0.0096)

Riquesa de coleòpters

# Filtrar les dades per coleòpters
  dades_coleoptera <- dades_finals[dades_finals$Ordre == "Coleoptera", ]
head(dades_coleoptera)  
## # A tibble: 6 × 26
##   Especimen Zona       Mes    `Mes num`   Any Trampa Tractament Identificacio   
##       <dbl> <chr>      <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>           
## 1        29 Campus     Juliol         7  2021      4 NS         <NA>            
## 2       158 Campus     Juliol         7  2021     10 NS         <NA>            
## 3       183 Campus     Juliol         7  2021      5 S          <NA>            
## 4       101 Franqueses Maig           5  2021      2 NS         <NA>            
## 5       150 Franqueses Maig           5  2021      5 NS         Chiqui bupresti…
## 6       151 Franqueses Maig           5  2021      5 NS         Chiqui bupresti…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
riquesa_x_trampa <- dades_coleoptera %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop") 

  proves_t_coleoptera <- riquesa_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Riquesa[Tractament == "S"]),
      mean_NS = mean(Riquesa[Tractament == "NS"]),
      t_test_result = tryCatch(
      t.test(Riquesa ~ Tractament)$p.value,
      error = function(e) NA
    )
      ,.groups = "drop")
  
  proves_t_coleoptera
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     1       1           NA    
## 2 Campus-Jun-2021        2.88    3.7          0.357
## 3 Campus-juliol-2020     1.33    1.25         0.851
## 4 Campus-maig-2020     NaN       2.44        NA    
## 5 Franqueses-Maig-2021   1.5     1            0.500
## 6 Sabadell-Juny-2021     1       1.33         0.423
  ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Riquesa x Zona i Tractaments",
       x = "Mostreig",
       y = "Riquesa")

ERROR: NOT ENOUGH OBSERVATIONS. ALGUNS MOSTREJOS NO TENEN COLEÒPTERS.

Riquesa de sírfids

# Filtrar les dades per sírfids
dades_sirfids <-dades_finals %>% filter(Familia=="Syrphidae")
head(dades_sirfids)  
## # A tibble: 6 × 26
##   Especimen Zona       Mes    `Mes num`   Any Trampa Tractament Identificacio   
##       <dbl> <chr>      <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>           
## 1       226 Campus     Juliol         7  2021      7 S          <NA>            
## 2       265 Campus     Juliol         7  2021     14 NS         <NA>            
## 3       277 Franqueses Maig           5  2021      4 NS         <NA>            
## 4       278 Franqueses Maig           5  2021      4 NS         <NA>            
## 5       279 Franqueses Maig           5  2021      4 NS         <NA>            
## 6       280 Franqueses Maig           5  2021      4 NS         Episyrphus balt…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
riquesa_x_trampa <- dades_sirfids %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop") 
  
  proves_t_sirfids <- riquesa_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Riquesa[Tractament == "S"]),
      mean_NS = mean(Riquesa[Tractament == "NS"]),
      t_test_result = tryCatch(
      t.test(Riquesa ~ Tractament)$p.value,
      error = function(e) NA
    )
      ,.groups = "drop")
  
  proves_t_sirfids
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl> <lgl>        
## 1 Campus-Juliol-2021        1    1    NA           
## 2 Campus-Jun-2021           1    1    NA           
## 3 Campus-juliol-2020        3    1    NA           
## 4 Campus-maig-2020        NaN    1.18 NA           
## 5 Franqueses-Maig-2021      1    3    NA           
## 6 Sabadell-Juny-2021      NaN    1    NA
ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Riquesa x Zona i Tractaments",
       x = "Mostreig",
       y = "Riquesa")

Hi ha poques dades per obtenir resultats del t-test.

Comparació de l’abundància per mostrejos i tractaments

Abundància total

abundancia_x_trampa <- dades_finals %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Abundancia = n(),.groups = "drop") 
  
  proves_t_abtotal <- abundancia_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Abundancia[Tractament == "S"]),
      mean_NS = mean(Abundancia[Tractament == "NS"]),
      t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(Abundancia ~ Tractament)$p.value,
                             NA)
      ,.groups = "drop")
  
  proves_t_abtotal
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021    14.1     4.5          0.130
## 2 Campus-Jun-2021       11.6    18.5          0.113
## 3 Campus-juliol-2020     8.89    8.54         0.904
## 4 Campus-maig-2020     NaN      15.4         NA    
## 5 Franqueses-Maig-2021   4.67    9            0.178
## 6 Sabadell-Juny-2021     7.29    5.8          0.606
 ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Abundancia x Zona i Tractaments",
       x = "Mostreig",
       y = "Abundancia") +
  scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
  scale_color_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
  theme_minimal()

No hi ha cap resultat significatiu (p-valor<0.05). No hi ha un patró clar de com varia l’abundància en funció del tractament.

Abundància d’abelles

# Filtrar les dades d'abelles
  dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
head(dades_abelles)  
## # A tibble: 6 × 26
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
abundancia_x_trampa <- dades_abelles %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Abundancia = n(),.groups = "drop") 
  
  proves_t_ababelles <- abundancia_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Abundancia[Tractament == "S"]),
      mean_NS = mean(Abundancia[Tractament == "NS"]),
      t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(Abundancia ~ Tractament)$p.value,
                             NA)
      ,.groups = "drop")
  
  proves_t_ababelles
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021    13.7     3.5         0.108 
## 2 Campus-Jun-2021        7.6     8.9         0.589 
## 3 Campus-juliol-2020     7.44    7.08        0.896 
## 4 Campus-maig-2020     NaN       8.68       NA     
## 5 Franqueses-Maig-2021   2       4.4         0.0902
## 6 Sabadell-Juny-2021     7       4.8         0.454
 ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Abundancia x Zona i Tractaments",
       x = "Mostreig",
       y = "Abundancia") +
  scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
  scale_color_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
  theme_minimal()

A Franqueses l’abundància d’abelles és significativament superior al tractament NS (p-valor=0.09015). En la resta de mostrejos el patró no és clar.

Abundància de coleòpters

# Filtrar les dades per coleòpters
  dades_coleoptera <- dades_finals[dades_finals$Ordre == "Coleoptera", ]
head(dades_coleoptera)  
## # A tibble: 6 × 26
##   Especimen Zona       Mes    `Mes num`   Any Trampa Tractament Identificacio   
##       <dbl> <chr>      <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>           
## 1        29 Campus     Juliol         7  2021      4 NS         <NA>            
## 2       158 Campus     Juliol         7  2021     10 NS         <NA>            
## 3       183 Campus     Juliol         7  2021      5 S          <NA>            
## 4       101 Franqueses Maig           5  2021      2 NS         <NA>            
## 5       150 Franqueses Maig           5  2021      5 NS         Chiqui bupresti…
## 6       151 Franqueses Maig           5  2021      5 NS         Chiqui bupresti…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
abundancia_x_trampa <- dades_coleoptera %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Abundancia = n(),.groups = "drop") 
  
  proves_t_abcoleoptera <- abundancia_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Abundancia[Tractament == "S"]),
      mean_NS = mean(Abundancia[Tractament == "NS"]),
      t_test_result = tryCatch(
      t.test(Abundancia ~ Tractament)$p.value,
      error = function(e) NA
    )
      ,.groups = "drop")
  
  proves_t_abcoleoptera
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     1       1           NA    
## 2 Campus-Jun-2021        4.12    9.2          0.175
## 3 Campus-juliol-2020     1.33    1.25         0.851
## 4 Campus-maig-2020     NaN       6.61        NA    
## 5 Franqueses-Maig-2021   2       4.67         0.496
## 6 Sabadell-Juny-2021     1       1.33         0.423
  ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Abundancia x Zona i Tractaments",
       x = "Mostreig",
       y = "Abundancia")

ERROR: NOT ENOUGH OBSERVATIONS. ALGUNS MOSTREJOS NO TENEN COLEÒPTERS.

Abundància de sírfids

# Filtrar les dades per sírfids
dades_sirfids <-dades_finals %>% filter(Familia=="Syrphidae")
head(dades_sirfids)  
## # A tibble: 6 × 26
##   Especimen Zona       Mes    `Mes num`   Any Trampa Tractament Identificacio   
##       <dbl> <chr>      <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>           
## 1       226 Campus     Juliol         7  2021      7 S          <NA>            
## 2       265 Campus     Juliol         7  2021     14 NS         <NA>            
## 3       277 Franqueses Maig           5  2021      4 NS         <NA>            
## 4       278 Franqueses Maig           5  2021      4 NS         <NA>            
## 5       279 Franqueses Maig           5  2021      4 NS         <NA>            
## 6       280 Franqueses Maig           5  2021      4 NS         Episyrphus balt…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
abundancia_x_trampa <- dades_sirfids %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Abundancia = n_distinct(ID),.groups = "drop") 
  
  proves_t_absirfids <- abundancia_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Abundancia[Tractament == "S"]),
      mean_NS = mean(Abundancia[Tractament == "NS"]),
      t_test_result = tryCatch(
      t.test(Abundancia ~ Tractament)$p.value,
      error = function(e) NA
    )
      ,.groups = "drop")
  
  proves_t_absirfids
## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl> <lgl>        
## 1 Campus-Juliol-2021        1    1    NA           
## 2 Campus-Jun-2021           1    1    NA           
## 3 Campus-juliol-2020        3    1    NA           
## 4 Campus-maig-2020        NaN    1.18 NA           
## 5 Franqueses-Maig-2021      1    3    NA           
## 6 Sabadell-Juny-2021      NaN    1    NA
  ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Abundancia x Zona i Tractaments",
       x = "Mostreig",
       y = "Abundancia")

No hi ha prous dades per obtenir bons resultats de la prova t-student.

DIAGRAMES D’EULER:

Diagrames d’Euler totals per cada zona

# Preparació de les dades: càlcul de l'abundància de cada espècie (ID) a cada mostreig i tractament (NS, S i intersecció)
taula_aux_euler <- as.data.frame.matrix(table(interaction(dades_finals$Zona,dades_finals$ID),dades_finals$Tractament)) %>% mutate(`NS&S`=pmin(NS,S)) 

head(taula_aux_euler)
##                                  NS S NS&S
## Campus.Acmaeodera cylindrica      5 1    1
## Franqueses.Acmaeodera cylindrica  0 0    0
## Sabadell.Acmaeodera cylindrica    0 0    0
## Campus.Amegilla sp.               1 0    0
## Franqueses.Amegilla sp.           0 0    0
## Sabadell.Amegilla sp.             0 0    0
# Reorganitzar la taula auxiliar perquè mostri el nom de les espècies (ID) trobades en cada mostreig i tractament
taula_euler_especies <- taula_aux_euler %>% mutate(Especie = sub("^[^.]+\\.", "", rownames(.)), Zona = sub("\\..*","",rownames(taula_aux_euler)))
taula_euler_especies <- taula_euler_especies %>% group_by(Zona) %>% summarise(  S = list(unique(Especie[ S > 0 ])),
    NS = list(unique(Especie[ NS > 0 ])),
   `NS&S` = list(unique(Especie[ `NS&S` > 0 ]))  )

head(taula_euler_especies)
## # A tibble: 3 × 4
##   Zona       S          NS         `NS&S`    
##   <chr>      <list>     <list>     <list>    
## 1 Campus     <chr [59]> <chr [98]> <chr [39]>
## 2 Franqueses <chr [11]> <chr [20]> <chr [5]> 
## 3 Sabadell   <chr [12]> <chr [16]> <chr [5]>
# Descarregar les llistes d'espècies en un excel
write.xlsx(taula_euler_especies, file.path(directorio, "2.0_Euler_Especies.xlsx"))
# Agrupació per mostrejos (Zona, Mes, Any)
euler_data <- taula_euler_especies %>%
  mutate(NS = map_dbl(NS, length)-map_dbl(`NS&S`, length),
         S = map_dbl(S, length)-map_dbl(`NS&S`, length),
         `NS&S` = map_dbl(`NS&S`, length))

head(euler_data)
## # A tibble: 3 × 4
##   Zona           S    NS `NS&S`
##   <chr>      <dbl> <dbl>  <dbl>
## 1 Campus        20    59     39
## 2 Franqueses     6    15      5
## 3 Sabadell       7    11      5
# Creació del diagrama d'Euler
apply_euler <- function(row){
  euler(c("NS"=row$NS,"S"=row$S,"NS&S"=row$`NS&S`))
}

euler_data <- euler_data %>% rowwise() %>% mutate(euler_result= list(apply_euler(cur_data())))
# Plot del diagrama d'Euler
quantities_font <- list(color = "white", size = 1.8, style = 2)

generate_euler_plots <- function(euler_result, zona, colors, label_size = 1.5){
  plot(euler_result, main = paste(zona), quantities = list(type = c("counts"), 
                         col = quantities_font$color, 
                         cex = quantities_font$size, 
                         font = quantities_font$style), fill = colors, 
       labels = list(col = "white", cex = label_size))
}

# Example of specifying a color palette
colors_palette <- c("#81b29a", "#e07a5f", "#3D405B")

# Use pmap to apply function with color argument
pmap(list(euler_data$euler_result, euler_data$Zona, list(colors_palette)), generate_euler_plots)
## [[1]]

## 
## [[2]]

## 
## [[3]]

Diagrames d’Euler d’abelles (amb dades_abelles) per cada zona

dades_abelles <- dades_abelles[!(dades_abelles$Mes == "maig" & dades_abelles$Any == 2020 & dades_abelles$Zona == "Campus"), ]
# Preparació de les dades: càlcul de l'abundància de cada espècie (ID) a cada mostreig i tractament (NS, S i intersecció)
taula_aux_eulerabelles <- as.data.frame.matrix(table(interaction(dades_abelles$Zona,dades_abelles$ID),dades_abelles$Tractament)) %>% mutate(`NS&S`=pmin(NS,S)) 

head(taula_aux_eulerabelles)
##                            NS S NS&S
## Campus.Amegilla sp.         1 0    0
## Franqueses.Amegilla sp.     0 0    0
## Sabadell.Amegilla sp.       0 0    0
## Campus.Andrena cinerea      0 0    0
## Franqueses.Andrena cinerea  1 0    0
## Sabadell.Andrena cinerea    0 0    0
# Reorganitzar la taula auxiliar perquè mostri el nom de les espècies (ID) trobades en cada mostreig i tractament
taula_euler_abelles <- taula_aux_eulerabelles %>% mutate(Especie = sub("^[^.]+\\.", "", rownames(.)), Zona = sub("\\..*","",rownames(taula_aux_eulerabelles)))
taula_euler_abelles <- taula_euler_abelles %>% group_by(Zona) %>% summarise(  S = list(unique(Especie[ S > 0 ])),
    NS = list(unique(Especie[ NS > 0 ])),
   `NS&S` = list(unique(Especie[ `NS&S` > 0 ]))  )

head(taula_euler_abelles)
## # A tibble: 3 × 4
##   Zona       S          NS         `NS&S`    
##   <chr>      <list>     <list>     <list>    
## 1 Campus     <chr [35]> <chr [47]> <chr [21]>
## 2 Franqueses <chr [5]>  <chr [12]> <chr [3]> 
## 3 Sabadell   <chr [10]> <chr [12]> <chr [5]>
# Descarregar les llistes d'espècies en un excel
write.xlsx(taula_euler_abelles, file.path(directorio, "Euler_Abelles.xlsx"))
# Agrupació per mostrejos (Zona, Mes, Any)
euler_data_abelles <- taula_euler_abelles %>%
  mutate(NS = map_dbl(NS, length)-map_dbl(`NS&S`, length),
         S = map_dbl(S, length)-map_dbl(`NS&S`, length),
         `NS&S` = map_dbl(`NS&S`, length))

head(euler_data_abelles)
## # A tibble: 3 × 4
##   Zona           S    NS `NS&S`
##   <chr>      <dbl> <dbl>  <dbl>
## 1 Campus        14    26     21
## 2 Franqueses     2     9      3
## 3 Sabadell       5     7      5
# Creació del diagrama d'Euler
apply_euler <- function(row){
  euler(c("NS"=row$NS,"S"=row$S,"NS&S"=row$`NS&S`))
}

euler_data_abelles <- euler_data_abelles %>% rowwise() %>% mutate(euler_result= list(apply_euler(cur_data())))
# Plot del diagrama d'Euler
quantities_font <- list(color = "white", size = 1.8, style = 2)
generate_euler_plots <- function(euler_result, zona, colors, label_size = 1.5){
  plot(euler_result, main = paste(zona), quantities = list(type = c("counts"), 
                         col = quantities_font$color, 
                         cex = quantities_font$size, 
                         font = quantities_font$style), fill = colors, 
       labels = list(col = "white", cex = label_size))
}

# Example of specifying a color palette
colors_palette <- c("#81b29a", "#e07a5f", "#3D405B")

# Use pmap to apply function with color argument
pmap(list(euler_data_abelles$euler_result, euler_data$Zona, list(colors_palette)), generate_euler_plots)
## [[1]]

## 
## [[2]]

## 
## [[3]]

ANÀLISIS DE COMPONENTS PRINCIPALS (PCA)

PCA total

# Preparar dades: taula de dades Shannon per mostreig i tractament (MT). Normalitzar les dades: nombre d'individus entre el nombre de pantraps del mostreig 
pca_data <- data_shannon_MT

# Calcular el nombre de trampes per cada mostreig i tractament
n_trampes <- dades_finals %>% group_by(Mostreig, Tractament) %>% summarise(numtrampes=n_distinct(Trampa)) %>% ungroup() %>%  mutate(row=paste(Mostreig,Tractament,sep="_")) %>% dplyr::select(-Mostreig,-Tractament)
## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.
# Unir les dues taules relacionant les files del mateix mostreig i tractament
pca_data$row <- rownames(pca_data)
pca_data <- pca_data  %>% inner_join(n_trampes, by=join_by(row))

# Dividir les abundàncies de cada ID entre el nombre de trampes de cada mostreig i tractament
num_columns <-ncol(pca_data) 
pca_data[, 1:(num_columns - 2)] <- pca_data[, 1:(num_columns - 2)] / pca_data[, num_columns]

# Eliminar l'última columna de la taula (numtrampes), ja no cal per fer la PCA
pca_data <- pca_data %>% dplyr::select(-numtrampes)

pca_data<- pca_data %>%
  mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2])
# Fer la PCA
pca_comps <- prcomp(pca_data %>% dplyr::select(-row,-Tractament))
# Mostrar la PCA
ggbiplot(pca_comps, labels = pca_data$row, groups = pca_data$Tractament, 
         ellipse = TRUE, var.axes = FALSE) +
  theme_minimal() +
  theme(text = element_text(size = 20))

PCA abelles

# Preparar dades: taula de dades Shannon d'abelles. Normalitzar les dades: nombre d'individus entre el nombre de pantraps del mostreig 
pca_data_abelles <- data_shannon_abelles

# Calcular el nombre de trampes per cada mostreig i tractament
n_trampes_abelles <- dades_finals %>%
  filter(Ordre == 'HymenopteraAbella') %>%
  group_by(Mostreig, Tractament) %>%
  summarise(numtrampes = n_distinct(Trampa)) %>%
  ungroup() %>%
  mutate(row = paste(Mostreig, Tractament, sep = "_")) %>%
  dplyr::select(-Mostreig, -Tractament)
## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.
# Unir les dues taules relacionant les files del mateix mostreig i tractament
pca_data_abelles$row <- rownames(pca_data_abelles)
pca_data_abelles <- pca_data_abelles  %>% inner_join(n_trampes_abelles, by=join_by(row))

# Dividir les abundàncies de cada ID entre el nombre de trampes de cada mostreig i tractament
num_columnsabe <-ncol(pca_data_abelles) 
pca_data_abelles[, 1:(num_columnsabe - 2)] <- pca_data_abelles[, 1:(num_columnsabe - 2)] / pca_data_abelles[, num_columnsabe]

# Eliminar l'última columna de la taula (numtrampes), ja no cal per fer la PCA
pca_data_abelles <- pca_data_abelles %>% dplyr::select(-numtrampes)

pca_data_abelles<- pca_data_abelles %>%
  mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2])
# Fer la PCA
pca_abelles <- prcomp(pca_data_abelles %>% dplyr::select(-row,-Tractament))
summary(pca_abelles)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     1.9981 1.1729 0.61490 0.52897 0.48323 0.38196 0.31303
## Proportion of Variance 0.6025 0.2076 0.05706 0.04223 0.03524 0.02202 0.01479
## Cumulative Proportion  0.6025 0.8101 0.86717 0.90940 0.94464 0.96665 0.98144
##                           PC8     PC9    PC10      PC11
## Standard deviation     0.2482 0.19318 0.15503 2.126e-16
## Proportion of Variance 0.0093 0.00563 0.00363 0.000e+00
## Cumulative Proportion  0.9907 0.99637 1.00000 1.000e+00
# Mostrar la PCA
ggbiplot(pca_abelles, labels = pca_data_abelles$row, groups = pca_data_abelles$Tractament, 
         ellipse = TRUE, var.axes = FALSE) +
  theme_minimal() +
  theme(text = element_text(size = 20))

PCA abelles Campus (per trampes. No sé ben bé què he fet)

# Preparar dades: taula de dades Shannon d'abelles per mostreig, tractament i trampa (dades NO normalitzades)
pca_data_campus <- data_shannon_MTTR %>% 
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  filter(Zona == 'Campus') %>% 
  mutate(Tractament = str_extract(rownames(.), "(?<=_)[^_]+(?=_)"))

#View(pca_data_campus)

# Fer la PCA
pca_campus <- prcomp(pca_data_campus %>% dplyr::select(-Zona,-Tractament))

# Mostrar la PCA
ggbiplot(pca_campus, labels = pca_data_campus$row, groups = pca_data_campus$Tractament, 
         ellipse = TRUE, var.axes = FALSE) +
  theme_minimal() +
  theme(text = element_text(size = 20))

NMDS, ANÀLISIS DE SIMILITUD (ANOSIM), ADONIS, I ESPÈCIES INDICADORES

NMDS pol·linitzadors Campus

dataNMDS_polcampus <- pca_data_campus %>% dplyr::select(-Zona, -Tractament)

# Extreure de la taula el tractament de cada mostreig: una nova columna
Metadata1 <- dataNMDS_polcampus %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Metadata1 <- Metadata1 %>% mutate(color=ifelse(Tractament=="NS","black","black"))
# NMDS Pol·linitzadors Campus: 
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)

set.seed(7)

# k=2: nombre de dimensions
NMDS_polcampus <- metaMDS(dataNMDS_polcampus, k=2)
## Wisconsin double standardization
## Run 0 stress 0.1731463 
## Run 1 stress 0.1733912 
## ... Procrustes: rmse 0.0970393  max resid 0.5241494 
## Run 2 stress 0.1727394 
## ... New best solution
## ... Procrustes: rmse 0.08088437  max resid 0.4655322 
## Run 3 stress 0.1732186 
## ... Procrustes: rmse 0.08060157  max resid 0.4245035 
## Run 4 stress 0.173302 
## Run 5 stress 0.1744142 
## Run 6 stress 0.1742466 
## Run 7 stress 0.1744684 
## Run 8 stress 0.1740094 
## Run 9 stress 0.1740286 
## Run 10 stress 0.1740606 
## Run 11 stress 0.1750397 
## Run 12 stress 0.1743163 
## Run 13 stress 0.1747583 
## Run 14 stress 0.1748468 
## Run 15 stress 0.1731277 
## ... Procrustes: rmse 0.01610656  max resid 0.0950575 
## Run 16 stress 0.1732682 
## Run 17 stress 0.1754394 
## Run 18 stress 0.1734799 
## Run 19 stress 0.1737735 
## Run 20 stress 0.1729431 
## ... Procrustes: rmse 0.05032583  max resid 0.2448404 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     16: no. of iterations >= maxit
##      4: stress ratio > sratmax
colors <- c("#74AA90", "#DB6443")

# Plot NMDS pol·linitzadors Campus
ordiplot(NMDS_polcampus,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata1$Tractament)
for (i in seq_along(unique_groups)) {
  ordihull(NMDS_polcampus, groups = Metadata1$Tractament, 
           show.groups = unique_groups[i], 
           draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_polcampus,display="species",col="black",air=1, cex = 1.2)

# Save the last plot
g <- recordPlot()

# Replay and save as PNG
width <- 17 # Width in inches
height <- 14  # Height in inches
res <- 500   # Resolution in dpi

# Replay and save as PNG with higher resolution and larger size
png("PolCampus.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()
## quartz_off_screen 
##                 2

Anàlisi de similitud (ANOSIM) Total

ano = anosim(dataNMDS_polcampus, Metadata1$Tractament, distance = "bray", permutations = 9999)
ano
## 
## Call:
## anosim(x = dataNMDS_polcampus, grouping = Metadata1$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.0432 
##       Significance: 0.1186 
## 
## Permutation: free
## Number of permutations: 9999

No signiticatiu (p valor > 0.05): la composició de les comunitats no és significativament diferent entre els tractaments S i NS. També s’observa al plot de NMDS: els polígons que corresponen als dos tractaments estan superposats. Això probablement sigui per haver fet la prova amb totes les zones de mostreig alhora.

ADONIS total

# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_Polcampus <- dataNMDS_polcampus %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

NMDStot.dist <- vegdist(wisconsin(sqrt(dataNMDS_polcampus)), k=2)

adonis(formula = NMDStot.dist~Tractament, data = Tractament_metadata_Polcampus, permutations = 999)$aov.tab
## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## Tractament  1     0.837 0.83702  1.9865 0.02453  0.004 **
## Residuals  79    33.286 0.42135         0.97547          
## Total      80    34.123                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha diferències significatives entre tractaments (p-valor = 0.25874), però la interacció entre zona i tractament explica un 45% de la variabilitat, i és un resultat significatiu (p-valor=0.03996).

Espècies indicadores dels tractaments S i NS total

indicadores = multipatt(dataNMDS_polcampus, Metadata1$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 137
##  Selected number of species: 7 
##  Number of species associated to 1 group: 7 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  4 
##                               stat p.value  
## Lasioglossum cf. pauxillum   0.253  0.0356 *
## Rhodanthidium septemdentatum 0.249  0.0181 *
## Panurgus dentipes            0.235  0.0413 *
## Oedemera nobilis             0.217  0.0305 *
## 
##  Group S  #sps.  3 
##                             stat p.value    
## Lasioglossum politum       0.339  0.0003 ***
## Lasioglossum glabriusculum 0.246  0.0124 *  
## Lasioglossum griseolum     0.198  0.0408 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

NMDS abelles Campus

dataNMDS_abecampus <- abelles_shannon_MTTR %>% 
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  filter(Zona == "Campus") %>% dplyr::select(-Zona)

# Extreure de la taula el tractament de cada mostreig: una nova columna
Metadata2 <- dataNMDS_abecampus %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Metadata2 <- Metadata2 %>% mutate(color=ifelse(Tractament=="NS","black","black"))
# NMDS Abelles Campus: 
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)

set.seed(7)

# k=2: nombre de dimensions
NMDS_abecampus <- metaMDS(dataNMDS_abecampus, k=2)
## Wisconsin double standardization
## Run 0 stress 0.1265607 
## Run 1 stress 0.1259806 
## ... New best solution
## ... Procrustes: rmse 0.04368297  max resid 0.2210331 
## Run 2 stress 0.1252628 
## ... New best solution
## ... Procrustes: rmse 0.05143489  max resid 0.2737022 
## Run 3 stress 0.1253757 
## ... Procrustes: rmse 0.03273285  max resid 0.1846694 
## Run 4 stress 0.1245387 
## ... New best solution
## ... Procrustes: rmse 0.04125358  max resid 0.2741621 
## Run 5 stress 0.1252543 
## Run 6 stress 0.1252991 
## Run 7 stress 0.1251473 
## Run 8 stress 0.1252565 
## Run 9 stress 0.125563 
## Run 10 stress 0.1258315 
## Run 11 stress 0.1249854 
## ... Procrustes: rmse 0.03920178  max resid 0.2852872 
## Run 12 stress 0.125473 
## Run 13 stress 0.1247876 
## ... Procrustes: rmse 0.03278848  max resid 0.2426446 
## Run 14 stress 0.1255484 
## Run 15 stress 0.1256229 
## Run 16 stress 0.1266907 
## Run 17 stress 0.1244476 
## ... New best solution
## ... Procrustes: rmse 0.02375557  max resid 0.1366641 
## Run 18 stress 0.1251499 
## Run 19 stress 0.125273 
## Run 20 stress 0.1256626 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: no. of iterations >= maxit
colors <- c("#74AA90", "#DB6443")

# Plot NMDS Abelles Campus
ordiplot(NMDS_abecampus,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata2$Tractament)
for (i in seq_along(unique_groups)) {
  ordihull(NMDS_abecampus, groups = Metadata2$Tractament, 
           show.groups = unique_groups[i], 
           draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_abecampus,display="species",col="black",air=1, cex = 1.2)

# Save the last plot
g <- recordPlot()

# Replay and save as PNG
width <- 17 # Width in inches
height <- 14  # Height in inches
res <- 500   # Resolution in dpi

# Replay and save as PNG with higher resolution and larger size
png("AbeCampus.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()
## quartz_off_screen 
##                 2

Anàlisi de similitud (ANOSIM) Abelles Campus

ano = anosim(dataNMDS_abecampus, Metadata2$Tractament, distance = "bray", permutations = 9999)
ano
## 
## Call:
## anosim(x = dataNMDS_abecampus, grouping = Metadata2$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.029 
##       Significance: 0.1956 
## 
## Permutation: free
## Number of permutations: 9999

No signiticatiu (p valor = 0.1956): la composició de les comunitats no és significativament diferent entre els tractaments S i NS.

ADONIS Abelles Campus

# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_Abecampus <- dataNMDS_abecampus %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

NMDStot.dist <- vegdist(wisconsin(sqrt(dataNMDS_abecampus)), k=2)

adonis(formula = NMDStot.dist~Tractament, data = Tractament_metadata_Abecampus, permutations = 999)$aov.tab
## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## Tractament  1     0.880 0.87965  2.1685 0.02705  0.006 **
## Residuals  78    31.641 0.40566         0.97295          
## Total      79    32.521                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El factor Tractament explica una pequeña proporción de la variabilidad total (R2 = 0.02705). La diferencia en la composición de las comunidades entre los grupos definidos por Tractament es estadísticamente significativa (valor p = 0.006). El estadístico F.Model de 2.1685 sugiere que hay una diferencia entre los grupos mayor que la esperada por el azar. Por lo tanto, a pesar de que la proporción de la varianza explicada es pequeña, la diferencia en la composición entre los grupos de Tractament es significativa desde un punto de vista estadístico.

Espècies indicadores abelles Campus

indicadores2 = multipatt(dataNMDS_abecampus, Metadata2$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores2)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 80
##  Selected number of species: 6 
##  Number of species associated to 1 group: 6 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  3 
##                               stat p.value  
## Lasioglossum cf. pauxillum   0.256  0.0362 *
## Rhodanthidium septemdentatum 0.252  0.0146 *
## Panurgus dentipes            0.241  0.0377 *
## 
##  Group S  #sps.  3 
##                             stat p.value    
## Lasioglossum politum       0.338  0.0005 ***
## Lasioglossum glabriusculum 0.244  0.0162 *  
## Lasioglossum griseolum     0.197  0.0435 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sí que hi ha espècies indicadores!!!!

NMDS pol·linitzadors Franqueses

dataNMDS_polfranqueses <- data_shannon_MTTR %>% 
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  filter(Zona == "Franqueses") %>% dplyr::select(-Zona)


# Extreure de la taula el tractament de cada mostreig: una nova columna
Metadata3 <- dataNMDS_polfranqueses %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Metadata3 <- Metadata3 %>% mutate(color=ifelse(Tractament=="NS","black","black"))
# NMDS: 
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)

set.seed(7)

# k=2: nombre de dimensions
NMDS_polfranqueses <- metaMDS(dataNMDS_polfranqueses, k=2)
## Wisconsin double standardization
## Run 0 stress 9.499321e-05 
## Run 1 stress 1.701317e-05 
## ... New best solution
## ... Procrustes: rmse 0.1971238  max resid 0.285099 
## Run 2 stress 0.04788483 
## Run 3 stress 8.301434e-05 
## ... Procrustes: rmse 0.1625411  max resid 0.2817144 
## Run 4 stress 9.344796e-05 
## ... Procrustes: rmse 0.2418918  max resid 0.3678843 
## Run 5 stress 8.989242e-05 
## ... Procrustes: rmse 0.1362289  max resid 0.221904 
## Run 6 stress 7.231338e-05 
## ... Procrustes: rmse 0.2054776  max resid 0.4848038 
## Run 7 stress 0 
## ... New best solution
## ... Procrustes: rmse 0.1785951  max resid 0.3935142 
## Run 8 stress 0 
## ... Procrustes: rmse 0.2030769  max resid 0.2694547 
## Run 9 stress 0.0477509 
## Run 10 stress 9.984028e-05 
## ... Procrustes: rmse 0.1689318  max resid 0.2881185 
## Run 11 stress 9.886818e-05 
## ... Procrustes: rmse 0.2345039  max resid 0.4475393 
## Run 12 stress 8.805342e-05 
## ... Procrustes: rmse 0.2059719  max resid 0.3315155 
## Run 13 stress 9.00509e-05 
## ... Procrustes: rmse 0.1985079  max resid 0.2915023 
## Run 14 stress 0.06801037 
## Run 15 stress 0.04833749 
## Run 16 stress 8.45766e-05 
## ... Procrustes: rmse 0.1123823  max resid 0.1615996 
## Run 17 stress 0 
## ... Procrustes: rmse 0.2343765  max resid 0.4118304 
## Run 18 stress 0 
## ... Procrustes: rmse 0.2541536  max resid 0.4202579 
## Run 19 stress 9.431488e-05 
## ... Procrustes: rmse 0.1992038  max resid 0.3022094 
## Run 20 stress 0.04775453 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      4: no. of iterations >= maxit
##     15: stress < smin
##      1: stress ratio > sratmax
## Warning in metaMDS(dataNMDS_polfranqueses, k = 2): stress is (nearly) zero: you
## may have insufficient data
## Warning in postMDS(out$points, dis, plot = max(0, plot - 1), ...): skipping
## half-change scaling: too few points below threshold
colors <- c("#74AA90", "#DB6443")

# Plot NMDS pol·linitzadors Campus
ordiplot(NMDS_polfranqueses,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata3$Tractament)
for (i in seq_along(unique_groups)) {
  ordihull(NMDS_polfranqueses, groups = Metadata3$Tractament, 
           show.groups = unique_groups[i], 
           draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_polfranqueses,display="species",col="black",air=1, cex = 1.2)

# Save the last plot
g <- recordPlot()

# Replay and save as PNG
width <- 17 # Width in inches
height <- 14  # Height in inches
res <- 500   # Resolution in dpi

# Replay and save as PNG with higher resolution and larger size
png("PolFranqueses.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()
## quartz_off_screen 
##                 2

Anàlisi de similitud (ANOSIM) Pol·linitzadors Franqueses

ano = anosim(dataNMDS_polfranqueses, Metadata3$Tractament, distance = "bray", permutations = 9999)
ano
## 
## Call:
## anosim(x = dataNMDS_polfranqueses, grouping = Metadata3$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.005128 
##       Significance: 0.4773 
## 
## Permutation: free
## Number of permutations: 9999

No signiticatiu (p valor = 0.4773): la composició de les comunitats no és significativament diferent entre els tractaments S i NS.

ADONIS Pol·linitzadors Franqueses

# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_polfranqueses <- dataNMDS_polfranqueses %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

head(Tractament_metadata_polfranqueses)
##                           Tractament
## Franqueses-Maig-2021_NS_1         NS
## Franqueses-Maig-2021_NS_2         NS
## Franqueses-Maig-2021_NS_3         NS
## Franqueses-Maig-2021_NS_4         NS
## Franqueses-Maig-2021_NS_5         NS
## Franqueses-Maig-2021_S_1           S
NMDStot.dist <- vegdist(wisconsin(sqrt(dataNMDS_polfranqueses)), k=2)

adonis(formula = NMDStot.dist~Tractament, data = Tractament_metadata_polfranqueses, permutations = 9999)$aov.tab
## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").
## Permutation: free
## Number of permutations: 9999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Tractament  1   0.39194 0.39194 0.89501 0.12981 0.6102
## Residuals   6   2.62749 0.43791         0.87019       
## Total       7   3.01943                 1.00000

No hi ha diferències significatives entre tractaments (p-valor = 0.6102).

Espècies indicadores Pol·linitzadors Franqueses

indicadores3 = multipatt(dataNMDS_polfranqueses, Metadata3$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores3)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 137
##  Selected number of species: 0 
##  Number of species associated to 1 group: 0 
## 
##  List of species associated to each combination: 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha espècies indicadores a Franqueses

NMDS abelles Franqueses

dataNMDS_abefranqueses <- abelles_shannon_MTTR %>% 
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  filter(Zona == "Franqueses") %>% dplyr::select(-Zona)

# Extreure de la taula el tractament de cada mostreig: una nova columna
Metadata4 <- dataNMDS_abefranqueses %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Metadata4 <- Metadata4 %>% mutate(color=ifelse(Tractament=="NS","black","black"))
# NMDS abelles Franqueses: 
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)

set.seed(7)

# k=2: nombre de dimensions
NMDS_abefranqueses <- metaMDS(dataNMDS_abefranqueses, k=2)
## Run 0 stress 0 
## Run 1 stress 0 
## ... Procrustes: rmse 0.1665426  max resid 0.251658 
## Run 2 stress 0 
## ... Procrustes: rmse 0.09615897  max resid 0.1451541 
## Run 3 stress 0 
## ... Procrustes: rmse 0.134433  max resid 0.2182935 
## Run 4 stress 9.652733e-05 
## ... Procrustes: rmse 0.1826024  max resid 0.2600087 
## Run 5 stress 9.726088e-05 
## ... Procrustes: rmse 0.1430138  max resid 0.2252209 
## Run 6 stress 3.030201e-05 
## ... Procrustes: rmse 0.1575005  max resid 0.2675314 
## Run 7 stress 0 
## ... Procrustes: rmse 0.1895275  max resid 0.3082301 
## Run 8 stress 0 
## ... Procrustes: rmse 0.117114  max resid 0.1693173 
## Run 9 stress 0 
## ... Procrustes: rmse 0.1365075  max resid 0.1737913 
## Run 10 stress 0 
## ... Procrustes: rmse 0.09004667  max resid 0.1559889 
## Run 11 stress 0 
## ... Procrustes: rmse 0.1075005  max resid 0.1820846 
## Run 12 stress 5.509784e-05 
## ... Procrustes: rmse 0.1815644  max resid 0.2208323 
## Run 13 stress 9.48924e-05 
## ... Procrustes: rmse 0.2171857  max resid 0.4526091 
## Run 14 stress 5.958231e-05 
## ... Procrustes: rmse 0.1243509  max resid 0.1690826 
## Run 15 stress 9.850827e-05 
## ... Procrustes: rmse 0.2171258  max resid 0.452336 
## Run 16 stress 4.866928e-05 
## ... Procrustes: rmse 0.2540161  max resid 0.507386 
## Run 17 stress 0 
## ... Procrustes: rmse 0.2351303  max resid 0.4343606 
## Run 18 stress 0 
## ... Procrustes: rmse 0.1837085  max resid 0.2811085 
## Run 19 stress 0 
## ... Procrustes: rmse 0.184819  max resid 0.3921269 
## Run 20 stress 9.072194e-05 
## ... Procrustes: rmse 0.1452401  max resid 0.1910837 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress < smin
## Warning in metaMDS(dataNMDS_abefranqueses, k = 2): stress is (nearly) zero: you
## may have insufficient data
## Warning in postMDS(out$points, dis, plot = max(0, plot - 1), ...): skipping
## half-change scaling: too few points below threshold
colors <- c("#74AA90", "#DB6443")

# Plot NMDS pol·linitzadors Campus
ordiplot(NMDS_abefranqueses,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata4$Tractament)
for (i in seq_along(unique_groups)) {
  ordihull(NMDS_abefranqueses, groups = Metadata4$Tractament, 
           show.groups = unique_groups[i], 
           draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_abefranqueses,display="species",col="black",air=1, cex = 1.2)

# Save the last plot
g <- recordPlot()

# Replay and save as PNG
width <- 17 # Width in inches
height <- 14  # Height in inches
res <- 500   # Resolution in dpi

# Replay and save as PNG with higher resolution and larger size
png("AbeFranqueses.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()
## quartz_off_screen 
##                 2

Anàlisi de similitud (ANOSIM) Abelles Franqueses

ano = anosim(dataNMDS_abefranqueses, Metadata4$Tractament, distance = "bray", permutations = 9999)
ano
## 
## Call:
## anosim(x = dataNMDS_abefranqueses, grouping = Metadata4$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.1128 
##       Significance: 0.2478 
## 
## Permutation: free
## Number of permutations: 9999

No signiticatiu (p valor = 0.2478): la composició de les comunitats no és significativament diferent entre els tractaments S i NS.

ADONIS Abelles Franqueses

# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_abefranqueses <- dataNMDS_abefranqueses %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

head(Tractament_metadata_abefranqueses)
##                           Tractament
## Franqueses-Maig-2021_NS_1         NS
## Franqueses-Maig-2021_NS_2         NS
## Franqueses-Maig-2021_NS_3         NS
## Franqueses-Maig-2021_NS_4         NS
## Franqueses-Maig-2021_NS_5         NS
## Franqueses-Maig-2021_S_1           S
NMDStot.dist <- vegdist(wisconsin(sqrt(dataNMDS_abefranqueses)), k=2)

adonis(formula = NMDStot.dist~Tractament, data = Tractament_metadata_abefranqueses, permutations = 9999)$aov.tab
## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").
## Permutation: free
## Number of permutations: 9999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Tractament  1   0.55855 0.55855  1.3907 0.18817 0.2104
## Residuals   6   2.40972 0.40162         0.81183       
## Total       7   2.96828                 1.00000

No hi ha diferències significatives entre tractaments (p-valor = 0.2104).

Espècies indicadores Abelles Franqueses

indicadores4 = multipatt(dataNMDS_abefranqueses, Metadata4$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores4)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 80
##  Selected number of species: 0 
##  Number of species associated to 1 group: 0 
## 
##  List of species associated to each combination: 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha espècies indicadores d’abelles a Franqueses.

NMDS Pol·linitzadors Sabadell

dataNMDS_polsabadell <- data_shannon_MTTR %>% 
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  filter(Zona == "Sabadell") %>% dplyr::select(-Zona)

# Extreure de la taula el tractament de cada mostreig: una nova columna
Metadata5 <- dataNMDS_polsabadell %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Metadata5 <- Metadata5 %>% mutate(color=ifelse(Tractament=="NS","black","black"))
# NMDS: 
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)

set.seed(7)

# k=2: nombre de dimensions
NMDS_polsabadell <- metaMDS(dataNMDS_polsabadell, k=2)
## Wisconsin double standardization
## Run 0 stress 0.1692459 
## Run 1 stress 0.1684865 
## ... New best solution
## ... Procrustes: rmse 0.1078519  max resid 0.2690531 
## Run 2 stress 0.178387 
## Run 3 stress 0.1720622 
## Run 4 stress 0.1746865 
## Run 5 stress 0.1693409 
## Run 6 stress 0.1951087 
## Run 7 stress 0.1783596 
## Run 8 stress 0.1782982 
## Run 9 stress 0.1684865 
## ... Procrustes: rmse 1.220837e-05  max resid 2.176932e-05 
## ... Similar to previous best
## Run 10 stress 0.1684865 
## ... Procrustes: rmse 3.457522e-06  max resid 6.800015e-06 
## ... Similar to previous best
## Run 11 stress 0.1684865 
## ... New best solution
## ... Procrustes: rmse 2.680791e-06  max resid 5.15384e-06 
## ... Similar to previous best
## Run 12 stress 0.188975 
## Run 13 stress 0.1683491 
## ... New best solution
## ... Procrustes: rmse 0.05383036  max resid 0.1315997 
## Run 14 stress 0.1925565 
## Run 15 stress 0.1693409 
## Run 16 stress 0.1720621 
## Run 17 stress 0.1783863 
## Run 18 stress 0.1955375 
## Run 19 stress 0.1900103 
## Run 20 stress 0.1757075 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     17: stress ratio > sratmax
##      2: scale factor of the gradient < sfgrmin
colors <- c("#74AA90", "#DB6443")

# Plot NMDS pol·linitzadors Campus
ordiplot(NMDS_polsabadell,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata5$Tractament)
for (i in seq_along(unique_groups)) {
  ordihull(NMDS_polsabadell, groups = Metadata5$Tractament, 
           show.groups = unique_groups[i], 
           draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_polsabadell,display="species",col="black",air=1, cex = 1.2)

# Save the last plot
g <- recordPlot()

# Replay and save as PNG
width <- 17 # Width in inches
height <- 14  # Height in inches
res <- 500   # Resolution in dpi

# Replay and save as PNG with higher resolution and larger size
png("PolSabadell.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()
## quartz_off_screen 
##                 2

Anàlisi de similitud (ANOSIM) Pol·linitzadors Sabadell

ano = anosim(dataNMDS_polsabadell, Metadata5$Tractament, distance = "bray", permutations = 9999)
ano
## 
## Call:
## anosim(x = dataNMDS_polsabadell, grouping = Metadata5$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.1244 
##       Significance: 0.8597 
## 
## Permutation: free
## Number of permutations: 9999

No signiticatiu (p valor = 0.8597): la composició de les comunitats no és significativament diferent entre els tractaments S i NS.

ADONIS Pol·linitzadors Sabadell

# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_polsabadell <- dataNMDS_polsabadell %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

head(Tractament_metadata_polsabadell)
##                          Tractament
## Sabadell-Juny-2021_NS_1          NS
## Sabadell-Juny-2021_NS_4          NS
## Sabadell-Juny-2021_NS_7          NS
## Sabadell-Juny-2021_NS_8          NS
## Sabadell-Juny-2021_NS_11         NS
## Sabadell-Juny-2021_S_2            S
NMDStot.dist <- vegdist(wisconsin(sqrt(dataNMDS_polsabadell)), k=2)

adonis(formula = NMDStot.dist~Tractament, data = Tractament_metadata_polsabadell, permutations = 9999)$aov.tab
## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").
## Permutation: free
## Number of permutations: 9999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Tractament  1    0.3582 0.35823   1.026 0.09305 0.4316
## Residuals  10    3.4916 0.34916         0.90695       
## Total      11    3.8498                 1.00000

No hi ha diferències significatives entre tractaments (p-valor = 0.4316).

Espècies indicadores Pol·linitzadors Sabadell

indicadores5 = multipatt(dataNMDS_polsabadell, Metadata5$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores5)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 137
##  Selected number of species: 0 
##  Number of species associated to 1 group: 0 
## 
##  List of species associated to each combination: 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha espècies indicadores de pol·linitzadors a Sabadell

NMDS Abelles Sabadell

dataNMDS_abesabadell <- abelles_shannon_MTTR %>% 
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  filter(Zona == "Sabadell") %>% dplyr::select(-Zona)

# Extreure de la taula el tractament de cada mostreig: una nova columna
Metadata6 <- dataNMDS_abesabadell %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Metadata6 <- Metadata6 %>% mutate(color=ifelse(Tractament=="NS","black","black"))
# NMDS abelles Sabadell: 
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)

set.seed(7)

# k=2: nombre de dimensions
NMDS_abesabadell <- metaMDS(dataNMDS_abesabadell, k=2)
## Wisconsin double standardization
## Run 0 stress 0.1316041 
## Run 1 stress 0.1469069 
## Run 2 stress 0.1595895 
## Run 3 stress 0.1335987 
## Run 4 stress 0.1657059 
## Run 5 stress 0.1316041 
## ... Procrustes: rmse 2.349789e-06  max resid 4.067505e-06 
## ... Similar to previous best
## Run 6 stress 0.1316041 
## ... Procrustes: rmse 2.165757e-05  max resid 3.932637e-05 
## ... Similar to previous best
## Run 7 stress 0.1316041 
## ... Procrustes: rmse 4.392255e-06  max resid 8.304797e-06 
## ... Similar to previous best
## Run 8 stress 0.1336417 
## Run 9 stress 0.1335987 
## Run 10 stress 0.1316041 
## ... Procrustes: rmse 2.786497e-05  max resid 5.137867e-05 
## ... Similar to previous best
## Run 11 stress 0.1337683 
## Run 12 stress 0.1316041 
## ... Procrustes: rmse 9.613328e-06  max resid 1.585993e-05 
## ... Similar to previous best
## Run 13 stress 0.1316041 
## ... Procrustes: rmse 3.396958e-06  max resid 6.825139e-06 
## ... Similar to previous best
## Run 14 stress 0.1695293 
## Run 15 stress 0.1330299 
## Run 16 stress 0.1335987 
## Run 17 stress 0.1316041 
## ... Procrustes: rmse 2.094775e-06  max resid 3.591374e-06 
## ... Similar to previous best
## Run 18 stress 0.169319 
## Run 19 stress 0.1673448 
## Run 20 stress 0.1732531 
## *** Best solution repeated 7 times
colors <- c("#74AA90", "#DB6443")

# Plot NMDS pol·linitzadors Campus
ordiplot(NMDS_abesabadell,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata6$Tractament)
for (i in seq_along(unique_groups)) {
  ordihull(NMDS_abesabadell, groups = Metadata6$Tractament, 
           show.groups = unique_groups[i], 
           draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_abesabadell,display="species",col="black",air=1, cex = 1.2)

# Save the last plot
g <- recordPlot()

# Replay and save as PNG
width <- 17 # Width in inches
height <- 14  # Height in inches
res <- 500   # Resolution in dpi

# Replay and save as PNG with higher resolution and larger size
png("AbeSabadell.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()
## quartz_off_screen 
##                 2

Anàlisi de similitud (ANOSIM) Abelles Sabadell

ano = anosim(dataNMDS_abesabadell, Metadata6$Tractament, distance = "bray", permutations = 9999)
ano
## 
## Call:
## anosim(x = dataNMDS_abesabadell, grouping = Metadata6$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.1622 
##       Significance: 0.9447 
## 
## Permutation: free
## Number of permutations: 9999

No signiticatiu (p valor = 0.9447): la composició de les comunitats no és significativament diferent entre els tractaments S i NS.

ADONIS Abelles Sabadell

# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_abesabadell <- dataNMDS_abesabadell %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

head(Tractament_metadata_abesabadell)
##                          Tractament
## Sabadell-Juny-2021_NS_1          NS
## Sabadell-Juny-2021_NS_4          NS
## Sabadell-Juny-2021_NS_7          NS
## Sabadell-Juny-2021_NS_8          NS
## Sabadell-Juny-2021_NS_11         NS
## Sabadell-Juny-2021_S_2            S
NMDStot.dist <- vegdist(wisconsin(sqrt(dataNMDS_abesabadell)), k=2)

adonis(formula = NMDStot.dist~Tractament, data = Tractament_metadata_abesabadell, permutations = 9999)$aov.tab
## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").
## Permutation: free
## Number of permutations: 9999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Tractament  1    0.2856 0.28561 0.86626 0.07972 0.5582
## Residuals  10    3.2970 0.32970         0.92028       
## Total      11    3.5826                 1.00000

No hi ha diferències significatives entre tractaments (p-valor = 0.5582).

Espècies indicadores Pol·linitzadors Sabadell

indicadores6 = multipatt(dataNMDS_abesabadell, Metadata6$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores6)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 80
##  Selected number of species: 0 
##  Number of species associated to 1 group: 0 
## 
##  List of species associated to each combination: 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha espècies d’abelles indicadores a Sabadell.

NMDS amb presència-absència (no abundàncies) ABELLES CAMPUS

NMDS_data_abellescampus_presabs<-dataNMDS_abecampus
NMDS_data_abellescampus_presabs[dataNMDS_abecampus>0]=1
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)

set.seed(7)

# k=2: nombre de dimensions
NMDS_presabs <- metaMDS(NMDS_data_abellescampus_presabs, k=2)
## Run 0 stress 0.1220077 
## Run 1 stress 0.1218802 
## ... New best solution
## ... Procrustes: rmse 0.05344429  max resid 0.2642867 
## Run 2 stress 0.1203376 
## ... New best solution
## ... Procrustes: rmse 0.07231672  max resid 0.3645329 
## Run 3 stress 0.1208759 
## Run 4 stress 0.119132 
## ... New best solution
## ... Procrustes: rmse 0.05694998  max resid 0.3191073 
## Run 5 stress 0.1191307 
## ... New best solution
## ... Procrustes: rmse 0.0306343  max resid 0.1370904 
## Run 6 stress 0.120841 
## Run 7 stress 0.1208037 
## Run 8 stress 0.1206953 
## Run 9 stress 0.1207318 
## Run 10 stress 0.1208075 
## Run 11 stress 0.1203058 
## Run 12 stress 0.120813 
## Run 13 stress 0.1202204 
## Run 14 stress 0.1203477 
## Run 15 stress 0.1230263 
## Run 16 stress 0.1232023 
## Run 17 stress 0.1202727 
## Run 18 stress 0.1220156 
## Run 19 stress 0.1206846 
## Run 20 stress 0.1212361 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     18: no. of iterations >= maxit
##      2: stress ratio > sratmax
colors <- c("#74AA90", "#DB6443")
# Plot NMDS total
ordiplot(NMDS_presabs,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata2$Tractament)
for (i in seq_along(unique_groups)) {
  ordihull(NMDS_presabs, groups = Metadata2$Tractament, 
           show.groups = unique_groups[i], 
           draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_presabs,display="species",col="black",air=0.01, cex = 1.2)


# Save the last plot
g <- recordPlot()

# Replay and save as PNG
width <- 17 # Width in inches
height <- 14  # Height in inches
res <- 500   # Resolution in dpi

# Replay and save as PNG with higher resolution and larger size
png("Presenciabsencia.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()

Anàlisi de similitud ANOSIM presència-absència

ano = anosim(NMDS_data_abellescampus_presabs, Metadata2$Tractament, distance = "bray", permutations = 9999)
ano
## 
## Call:
## anosim(x = NMDS_data_abellescampus_presabs, grouping = Metadata2$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.01252 
##       Significance: 0.3341 
## 
## Permutation: free
## Number of permutations: 9999

Les composicions de les comunitats dels dos tractaments son molt similars, i el resultat no és significatiu (0.3341).

ADONIS presència-absència

# Extreure de la taula el tractament de cada mostreig: una nova columna
ADONIS_Tractament_metadata <- NMDS_data_abellescampus_presabs %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

ADONIS_Tractament_metadata
##                          Tractament
## Campus-juliol-2020_NS_2          NS
## Campus-juliol-2020_NS_3          NS
## Campus-juliol-2020_NS_4          NS
## Campus-juliol-2020_NS_6          NS
## Campus-juliol-2020_NS_8          NS
## Campus-juliol-2020_NS_10         NS
## Campus-juliol-2020_NS_11         NS
## Campus-juliol-2020_NS_12         NS
## Campus-juliol-2020_NS_13         NS
## Campus-juliol-2020_NS_14         NS
## Campus-juliol-2020_NS_15         NS
## Campus-juliol-2020_NS_16         NS
## Campus-juliol-2020_NS_20         NS
## Campus-juliol-2020_S_1            S
## Campus-juliol-2020_S_5            S
## Campus-juliol-2020_S_7            S
## Campus-juliol-2020_S_9            S
## Campus-juliol-2020_S_17           S
## Campus-juliol-2020_S_18           S
## Campus-juliol-2020_S_19           S
## Campus-juliol-2020_S_21           S
## Campus-juliol-2020_S_22           S
## Campus-Juliol-2021_NS_2          NS
## Campus-Juliol-2021_NS_4          NS
## Campus-Juliol-2021_NS_6          NS
## Campus-Juliol-2021_NS_8          NS
## Campus-Juliol-2021_NS_10         NS
## Campus-Juliol-2021_NS_12         NS
## Campus-Juliol-2021_NS_14         NS
## Campus-Juliol-2021_NS_16         NS
## Campus-Juliol-2021_NS_18         NS
## Campus-Juliol-2021_NS_20         NS
## Campus-Juliol-2021_S_1            S
## Campus-Juliol-2021_S_3            S
## Campus-Juliol-2021_S_5            S
## Campus-Juliol-2021_S_7            S
## Campus-Juliol-2021_S_9            S
## Campus-Juliol-2021_S_11           S
## Campus-Juliol-2021_S_13           S
## Campus-Juliol-2021_S_15           S
## Campus-Juliol-2021_S_17           S
## Campus-Jun-2021_NS_2             NS
## Campus-Jun-2021_NS_4             NS
## Campus-Jun-2021_NS_6             NS
## Campus-Jun-2021_NS_8             NS
## Campus-Jun-2021_NS_10            NS
## Campus-Jun-2021_NS_12            NS
## Campus-Jun-2021_NS_14            NS
## Campus-Jun-2021_NS_16            NS
## Campus-Jun-2021_NS_18            NS
## Campus-Jun-2021_NS_20            NS
## Campus-Jun-2021_S_1               S
## Campus-Jun-2021_S_3               S
## Campus-Jun-2021_S_5               S
## Campus-Jun-2021_S_7               S
## Campus-Jun-2021_S_9               S
## Campus-Jun-2021_S_11              S
## Campus-Jun-2021_S_13              S
## Campus-Jun-2021_S_15              S
## Campus-Jun-2021_S_17              S
## Campus-Jun-2021_S_19              S
## Campus-maig-2020_NS_1            NS
## Campus-maig-2020_NS_2            NS
## Campus-maig-2020_NS_3            NS
## Campus-maig-2020_NS_4            NS
## Campus-maig-2020_NS_5            NS
## Campus-maig-2020_NS_6            NS
## Campus-maig-2020_NS_7            NS
## Campus-maig-2020_NS_8            NS
## Campus-maig-2020_NS_9            NS
## Campus-maig-2020_NS_10           NS
## Campus-maig-2020_NS_11           NS
## Campus-maig-2020_NS_12           NS
## Campus-maig-2020_NS_13           NS
## Campus-maig-2020_NS_14           NS
## Campus-maig-2020_NS_15           NS
## Campus-maig-2020_NS_16           NS
## Campus-maig-2020_NS_17           NS
## Campus-maig-2020_NS_18           NS
## Campus-maig-2020_NS_19           NS
NMDSabe.dist <- vegdist(wisconsin(sqrt(NMDS_data_abellescampus_presabs)), k=2)

adonis(formula = NMDSabe.dist~Tractament, data = ADONIS_Tractament_metadata, permutations = 1000)$aov.tab
## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").
## Permutation: free
## Number of permutations: 1000
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2  Pr(>F)  
## Tractament  1     0.818 0.81797  2.1994 0.02742 0.01299 *
## Residuals  78    29.009 0.37191         0.97258          
## Total      79    29.827                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El tractament explica un 2,724% de la variabilitat entre les comunitats, i el resultat és significatiu (p-valor = 0.01299).

Espècies indicadores presència-absència

indicadores = multipatt(NMDS_data_abellescampus_presabs, Metadata2$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores, alpha = 0.05)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 80
##  Selected number of species: 3 
##  Number of species associated to 1 group: 3 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  2 
##                               stat p.value  
## Lasioglossum cf. pauxillum   0.289  0.0270 *
## Rhodanthidium septemdentatum 0.287  0.0275 *
## 
##  Group S  #sps.  1 
##                       stat p.value  
## Lasioglossum politum 0.312  0.0109 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hi ha espècies indicadores de cada tractament!

MODELS LINEALS MIXTOS D’ANÀLISIS DESCRIPTIVES:

Abundància de Lasioglossums

# Crear matriu amb només les columnes de Lasioglossums

noms_columnes <- colnames(metadata_shannon)
genere <- "Lasioglossum"

columnes_lasioglossum <- noms_columnes[grep(paste0("^", genere), noms_columnes)]

dades_lasioglossum <- metadata_shannon[, columnes_lasioglossum]

# Agrupar totes les columnes en una sola (abundància total de Lasioglossums en cada trampa)
dades_lasioglossum$Abundancia_Total <- rowSums(dades_lasioglossum)
dades_lasioglossum_ok <- data.frame(Abundancia_Total = dades_lasioglossum$Abundancia_Total)
rownames(dades_lasioglossum_ok) <- rownames(dades_lasioglossum)

# View(dades_lasioglossum_ok)

# Extreure columnes de zona i tractament
lasioglossums_zona_tractament <- dades_lasioglossum_ok %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  mutate(Tractament = str_extract(rownames(.), "(?<=_)[A-Z]+(?=_\\d)"))

# View(lasioglossums_zona_tractament)

Model lineal mixt (tractament com a variable explicativa, i zona com a variable aleatòria)

# install.packages("lme4")
# install.packages("Matrix")

model_lasioglossum <- lmer(Abundancia_Total ~ Tractament + (1|Zona), data = lasioglossums_zona_tractament)
summary(model_lasioglossum)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Abundancia_Total ~ Tractament + (1 | Zona)
##    Data: lasioglossums_zona_tractament
## 
## REML criterion at convergence: 642.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1440 -0.4778 -0.1408  0.1962  7.2804 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Zona     (Intercept)  1.557   1.248   
##  Residual             35.226   5.935   
## Number of obs: 101, groups:  Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)   
## (Intercept)    1.926      1.172  2.051   1.643  0.23908   
## TractamentS    3.954      1.226 98.906   3.224  0.00171 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.429
plot(model_lasioglossum)

Quan el tractament és NS, l’abundància de Lasioglossums és superior a 0 (Intercept: marginalment significatiu). Quan el tractament és S, s’espera que l’abundància de lasioglossums augmenti 3,954 respecte el tractament NS, i és un resultat molt significatiu (p-valor = 0.00171). Els Lasioglossums son més abundants a les zones segades (concorda amb algunes de les espècies “indicadores” de S!)

Model lineal (Tractament com a variable explicativa, sense considerar la zona)

# Convertir la variable Tractament a factor
lasioglossums_zona_tractament$Tractament <- factor(lasioglossums_zona_tractament$Tractament)

# Ajustar el model lineal
model <- lm(Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament)
summary(model)
## 
## Call:
## lm(formula = Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.421 -2.619 -1.421  1.381 43.579 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.619      0.754   3.473 0.000764 ***
## TractamentS    3.802      1.229   3.093 0.002576 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.985 on 99 degrees of freedom
## Multiple R-squared:  0.08811,    Adjusted R-squared:  0.0789 
## F-statistic: 9.566 on 1 and 99 DF,  p-value: 0.002576
# Crear el gráfico
plot(Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament, pch = 16)

Diversitat de Shannon d’abelles (zona com a variable aleatòria)

# Crear un dataframe de la diversitat de shannon d'abelles (calculada abans)
taula_shannon_abelles <- data.frame(Shannon = abelles_shannondiv)
rownames(taula_shannon_abelles) <- rownames(abelles_shannon_MTTR)

# Extreure dues noves columnes: tractament i zona
taula_shannon_abelles <- taula_shannon_abelles %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2])

print(taula_shannon_abelles)
##                             Shannon       Zona Tractament
## Campus-juliol-2020_NS_2   1.3862944     Campus         NS
## Campus-juliol-2020_NS_3   1.6674619     Campus         NS
## Campus-juliol-2020_NS_4   1.6417347     Campus         NS
## Campus-juliol-2020_NS_6   1.3296613     Campus         NS
## Campus-juliol-2020_NS_8   1.6674619     Campus         NS
## Campus-juliol-2020_NS_10  0.7963116     Campus         NS
## Campus-juliol-2020_NS_11  0.0000000     Campus         NS
## Campus-juliol-2020_NS_12  1.5607104     Campus         NS
## Campus-juliol-2020_NS_13  1.6434177     Campus         NS
## Campus-juliol-2020_NS_14  0.7549968     Campus         NS
## Campus-juliol-2020_NS_15  1.5498260     Campus         NS
## Campus-juliol-2020_NS_16  0.6931472     Campus         NS
## Campus-juliol-2020_NS_20  1.7478681     Campus         NS
## Campus-juliol-2020_S_1    0.6931472     Campus          S
## Campus-juliol-2020_S_5    1.3517840     Campus          S
## Campus-juliol-2020_S_7    0.5505734     Campus          S
## Campus-juliol-2020_S_9    0.0000000     Campus          S
## Campus-juliol-2020_S_17   1.0986123     Campus          S
## Campus-juliol-2020_S_18   1.3862944     Campus          S
## Campus-juliol-2020_S_19   0.6931472     Campus          S
## Campus-juliol-2020_S_21   1.1189689     Campus          S
## Campus-juliol-2020_S_22   0.6108643     Campus          S
## Campus-Juliol-2021_NS_2   0.0000000     Campus         NS
## Campus-Juliol-2021_NS_4   0.0000000     Campus         NS
## Campus-Juliol-2021_NS_6   0.5623351     Campus         NS
## Campus-Juliol-2021_NS_8   1.0042425     Campus         NS
## Campus-Juliol-2021_NS_10  0.0000000     Campus         NS
## Campus-Juliol-2021_NS_12  1.0549202     Campus         NS
## Campus-Juliol-2021_NS_14  1.0397208     Campus         NS
## Campus-Juliol-2021_NS_16  0.0000000     Campus         NS
## Campus-Juliol-2021_NS_18  1.3862944     Campus         NS
## Campus-Juliol-2021_NS_20  0.6931472     Campus         NS
## Campus-Juliol-2021_S_1    0.6931472     Campus          S
## Campus-Juliol-2021_S_3    0.6365142     Campus          S
## Campus-Juliol-2021_S_5    1.8662160     Campus          S
## Campus-Juliol-2021_S_7    1.4369665     Campus          S
## Campus-Juliol-2021_S_9    1.1973401     Campus          S
## Campus-Juliol-2021_S_11   1.4388830     Campus          S
## Campus-Juliol-2021_S_13   0.0000000     Campus          S
## Campus-Juliol-2021_S_15   0.6931472     Campus          S
## Campus-Juliol-2021_S_17   0.5091373     Campus          S
## Campus-Jun-2021_NS_2      1.4735024     Campus         NS
## Campus-Jun-2021_NS_4      1.0986123     Campus         NS
## Campus-Jun-2021_NS_6      0.8675632     Campus         NS
## Campus-Jun-2021_NS_8      2.3393717     Campus         NS
## Campus-Jun-2021_NS_10     1.3862944     Campus         NS
## Campus-Jun-2021_NS_12     2.0947290     Campus         NS
## Campus-Jun-2021_NS_14     1.4925477     Campus         NS
## Campus-Jun-2021_NS_16     1.2130076     Campus         NS
## Campus-Jun-2021_NS_18     1.5607104     Campus         NS
## Campus-Jun-2021_NS_20     0.0000000     Campus         NS
## Campus-Jun-2021_S_1       1.6094379     Campus          S
## Campus-Jun-2021_S_3       1.5498260     Campus          S
## Campus-Jun-2021_S_5       1.0397208     Campus          S
## Campus-Jun-2021_S_7       0.5623351     Campus          S
## Campus-Jun-2021_S_9       1.5810938     Campus          S
## Campus-Jun-2021_S_11      1.5595812     Campus          S
## Campus-Jun-2021_S_13      1.4388830     Campus          S
## Campus-Jun-2021_S_15      0.4101163     Campus          S
## Campus-Jun-2021_S_17      1.3321790     Campus          S
## Campus-Jun-2021_S_19      1.2424533     Campus          S
## Campus-maig-2020_NS_1     0.6682485     Campus         NS
## Campus-maig-2020_NS_2     0.9611277     Campus         NS
## Campus-maig-2020_NS_3     0.6365142     Campus         NS
## Campus-maig-2020_NS_4     1.6094379     Campus         NS
## Campus-maig-2020_NS_5     0.6365142     Campus         NS
## Campus-maig-2020_NS_6     1.8845210     Campus         NS
## Campus-maig-2020_NS_7     0.6931472     Campus         NS
## Campus-maig-2020_NS_8     1.3862944     Campus         NS
## Campus-maig-2020_NS_9     1.3321790     Campus         NS
## Campus-maig-2020_NS_10    0.9502705     Campus         NS
## Campus-maig-2020_NS_11    1.3862944     Campus         NS
## Campus-maig-2020_NS_12    0.6931472     Campus         NS
## Campus-maig-2020_NS_13    1.3114313     Campus         NS
## Campus-maig-2020_NS_14    1.3296613     Campus         NS
## Campus-maig-2020_NS_15    1.4978661     Campus         NS
## Campus-maig-2020_NS_16    0.9502705     Campus         NS
## Campus-maig-2020_NS_17    0.6365142     Campus         NS
## Campus-maig-2020_NS_18    1.3862944     Campus         NS
## Campus-maig-2020_NS_19    0.9251291     Campus         NS
## Franqueses-Maig-2021_NS_1 0.6931472 Franqueses         NS
## Franqueses-Maig-2021_NS_2 0.5004024 Franqueses         NS
## Franqueses-Maig-2021_NS_3 0.6931472 Franqueses         NS
## Franqueses-Maig-2021_NS_4 1.6674619 Franqueses         NS
## Franqueses-Maig-2021_NS_5 1.0986123 Franqueses         NS
## Franqueses-Maig-2021_S_1  0.0000000 Franqueses          S
## Franqueses-Maig-2021_S_2  1.0986123 Franqueses          S
## Franqueses-Maig-2021_S_4  0.6931472 Franqueses          S
## Sabadell-Juny-2021_NS_1   1.7478681   Sabadell         NS
## Sabadell-Juny-2021_NS_4   0.0000000   Sabadell         NS
## Sabadell-Juny-2021_NS_7   1.3862944   Sabadell         NS
## Sabadell-Juny-2021_NS_8   1.2424533   Sabadell         NS
## Sabadell-Juny-2021_NS_11  1.3862944   Sabadell         NS
## Sabadell-Juny-2021_S_2    0.6365142   Sabadell          S
## Sabadell-Juny-2021_S_3    0.5623351   Sabadell          S
## Sabadell-Juny-2021_S_5    1.4995094   Sabadell          S
## Sabadell-Juny-2021_S_6    0.6931472   Sabadell          S
## Sabadell-Juny-2021_S_9    0.6931472   Sabadell          S
## Sabadell-Juny-2021_S_10   0.5004024   Sabadell          S
## Sabadell-Juny-2021_S_12   1.0986123   Sabadell          S
# Fer el model
model_shannonabelles <- lmer(Shannon ~ Tractament + (1|Zona), data = taula_shannon_abelles)
## boundary (singular) fit: see help('isSingular')
summary(model_shannonabelles)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Tractament + (1 | Zona)
##    Data: taula_shannon_abelles
## 
## REML criterion at convergence: 162.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0268 -0.7267  0.1090  0.7752  2.3611 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  Zona     (Intercept) 6.132e-18 2.476e-09
##  Residual             2.842e-01 5.331e-01
## Number of obs: 100, groups:  Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  1.08059    0.06771 98.00000  15.959   <2e-16 ***
## TractamentS -0.13912    0.10984 98.00000  -1.267    0.208    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.616
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Quan el tractament és NS, el valor de l’índex de Shannon és significativament superior a 0 (Intercept). Quan el tractament és S, s’espera que l’índex de Shannon disminueixi 0,4605 respecte el tractament NS. És un resultat significatiu (p-valor = 0,011055). Per tant, la diversitat de Shannon en abelles és significativament inferior en el tractament segat.

# Plot del model
plot(fitted(model_shannonabelles), resid(model_shannonabelles), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Riquesa total

# Preparar les dades
riquesa_aux <- as.data.frame(riquesa)
rownames(riquesa_aux) <- paste(riquesa_aux$Mostreig, riquesa_aux$Tractament, sep = "_")
riquesa_aux <- subset(riquesa_aux, select = -Mostreig)

# Extreure columna de la zona
taula_riquesa <- riquesa_aux %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1])

print(taula_riquesa)
##                         Tractament riquesa       Zona
## Campus-Juliol-2021_NS           NS      20     Campus
## Campus-Juliol-2021_S             S      21     Campus
## Campus-Jun-2021_NS              NS      45     Campus
## Campus-Jun-2021_S                S      34     Campus
## Campus-juliol-2020_NS           NS      39     Campus
## Campus-juliol-2020_S             S      24     Campus
## Campus-maig-2020_NS             NS      45     Campus
## Franqueses-Maig-2021_NS         NS      20 Franqueses
## Franqueses-Maig-2021_S           S      11 Franqueses
## Sabadell-Juny-2021_NS           NS      16   Sabadell
## Sabadell-Juny-2021_S             S      12   Sabadell
# Fer el model
model_riquesa <- lmer(riquesa ~ Tractament + (1|Zona), data = taula_riquesa)
summary(model_riquesa)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: riquesa ~ Tractament + (1 | Zona)
##    Data: taula_riquesa
## 
## REML criterion at convergence: 70.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7954 -0.4064 -0.1774  0.6667  1.1012 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Zona     (Intercept) 86.50    9.301   
##  Residual             74.49    8.631   
## Number of obs: 11, groups:  Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)  
## (Intercept)   26.171      6.627  2.885   3.949   0.0311 *
## TractamentS   -9.501      5.237  7.267  -1.814   0.1109  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.373

TractamentS: quan el tractament és S la riquesa disminueix (-9.501), però no és un resultat significatiu.

# Plot del model
plot(fitted(model_riquesa), resid(model_riquesa), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Riquesa PER TRAMPA d’abelles al Campus (data de mostreig com a factor aleatori)

# Preparar les dades
dades_abellescampus <- dades_abelles[dades_abelles$Zona == "Campus", ]

riquesatrampa_abellescampus <- dades_abellescampus %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop")

riquesa_abellescampus <- as.data.frame(riquesatrampa_abellescampus)
rownames(riquesa_abellescampus) <- paste(riquesa_abellescampus$Mostreig, riquesa_abellescampus$Tractament, riquesa_abellescampus$Trampa, sep = "_")
riquesa_abellescampus <- subset(riquesa_abellescampus, select = -Mostreig)

# Extreure columna de la data de mostreig:
riquesa_abellescampus <- riquesa_abellescampus %>%
  mutate(Data_Mostreig = str_extract(rownames(.), "(?<=-)[^-_]+-[0-9]+"))

head(riquesa_abellescampus)
##                          Tractament Trampa Riquesa Data_Mostreig
## Campus-Juliol-2021_NS_2          NS      2       1   Juliol-2021
## Campus-Juliol-2021_NS_4          NS      4       1   Juliol-2021
## Campus-Juliol-2021_NS_6          NS      6       2   Juliol-2021
## Campus-Juliol-2021_NS_8          NS      8       3   Juliol-2021
## Campus-Juliol-2021_NS_10         NS     10       1   Juliol-2021
## Campus-Juliol-2021_NS_12         NS     12       3   Juliol-2021
# Fer el model
model_riquesa_abellescampus <- lmer(Riquesa ~ Tractament + (1|Data_Mostreig), data = riquesa_abellescampus)
summary(model_riquesa_abellescampus)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Riquesa ~ Tractament + (1 | Data_Mostreig)
##    Data: riquesa_abellescampus
## 
## REML criterion at convergence: 265
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8026 -0.7445 -0.1710  0.4959  2.9434 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  Data_Mostreig (Intercept) 0.6885   0.8297  
##  Residual                  4.4396   2.1070  
## Number of obs: 61, groups:  Data_Mostreig, 3
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)   
## (Intercept)   4.0378     0.6041  2.8500   6.684  0.00801 **
## TractamentS  -0.3864     0.5427 57.1316  -0.712  0.47932   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.414

En el tractament S la riquesa d’abelles per trampa disminueix, però no és significatiu.

# Plot del model
plot(fitted(model_riquesa_abellescampus), resid(model_riquesa_abellescampus), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Riquesa PER MOSTREIG d’abelles al Campus (data de mostreig com a factor aleatori)

# Preparar les dades
riquesa_aux <- dades_abellescampus %>%
  group_by(Mostreig, Tractament) %>%
  summarise(riquesa = n_distinct(ID))
## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.
riquesa_aux <- as.data.frame(riquesa_aux)
rownames(riquesa_aux) <- paste(riquesa_aux$Mostreig, riquesa_aux$Tractament, sep = "_")
riquesa_aux <- subset(riquesa_aux, select = -Mostreig)

# Extreure columna de la data de mostreig:
riquesa_aux <- riquesa_aux %>%
  mutate(Data_Mostreig = str_extract(rownames(.), "(?<=-)[^-_]+-[0-9]+"))

riquesa_aux
##                       Tractament riquesa Data_Mostreig
## Campus-Juliol-2021_NS         NS      12   Juliol-2021
## Campus-Juliol-2021_S           S      17   Juliol-2021
## Campus-Jun-2021_NS            NS      25      Jun-2021
## Campus-Jun-2021_S              S      19      Jun-2021
## Campus-juliol-2020_NS         NS      27   juliol-2020
## Campus-juliol-2020_S           S      14   juliol-2020
# Fer el model
model_riquesa_abellescampus2 <- lmer(riquesa ~ Tractament + (1|Data_Mostreig), data = riquesa_aux)
## boundary (singular) fit: see help('isSingular')
summary(model_riquesa_abellescampus2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: riquesa ~ Tractament + (1 | Data_Mostreig)
##    Data: riquesa_aux
## 
## REML criterion at convergence: 27.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5484 -0.3180  0.2212  0.5530  0.9401 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  Data_Mostreig (Intercept)  0.00    0.000   
##  Residual                  36.33    6.028   
## Number of obs: 6, groups:  Data_Mostreig, 3
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)   
## (Intercept)   21.333      3.480  4.000   6.130  0.00359 **
## TractamentS   -4.667      4.922  4.000  -0.948  0.39672   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.707
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

TractamentS: quan el tractament és S la riquesa disminueix (en 5.583), però NO és un resultat significatiu (p-valor = 0.24677).

# Plot del model
plot(fitted(model_riquesa_abellescampus2), resid(model_riquesa_abellescampus2), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Diversitat (Shannon) d’abelles al Campus (data de mostreig com a factor aleatori)

# Crear taula auxiliar: abundàncies d'abelles per ID, per cada mostreig del Campus

# Filtrar dades per Campus
dades_abellescampus <- dades_finals[dades_finals$Zona == "Campus", ]

# Crear taula de dades d'abundància només per abelles
abundanciabelles_ID2 <- dades_abellescampus %>% 
  mutate(Mostreig = paste(Zona, Mes, Any, sep = "-")) %>% 
  filter(Ordre == 'HymenopteraAbella') %>% 
  group_by(Mostreig, Tractament, ID) %>%
  summarise(Numero_total = n())
## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.
# View(abundanciabelles_ID)

#Crear una taula auxiliar de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux2 <- cast(abundanciabelles_ID2, Mostreig + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux2 <- as.data.frame(taula_aux2)
taula_aux2[is.na(taula_aux2)] <- 0

#Agrupar Mostreig i tractament en una nova variable (Agrupacio)
data_shannon_abellescampus <- taula_aux2 %>% 
  mutate(Agrupacio = paste(Mostreig, Tractament, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament)

#Definir els rownames amb la variable Agrupacio (mostreig i tractament)
rownames(data_shannon_abellescampus) <- data_shannon_abellescampus$Agrupacio
data_shannon_abellescampus <- data_shannon_abellescampus %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

# head(data_shannon_abellescampus)
# Càlcul de l'índex de Shannon i preparació del dataframe
shannondiv_abellescampus <- diversity(data_shannon_abellescampus)

# Crear un dataframe de la diversitat de shannon d'abelles
taula_shannon_abellescampus <- data.frame(Shannon = shannondiv_abellescampus)
rownames(taula_shannon_abellescampus) <- rownames(data_shannon_abellescampus)

# Extreure dues noves columnes: tractament i data mostreig
taula_shannon_abellescampus <- taula_shannon_abellescampus %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% 
  mutate(Data_Mostreig = str_extract(rownames(.), "(?<=-)[^-_]+-[0-9]+"))

print(taula_shannon_abellescampus)
##                        Shannon Tractament Data_Mostreig
## Campus-juliol-2020_NS 2.611229         NS   juliol-2020
## Campus-juliol-2020_S  1.905899          S   juliol-2020
## Campus-Juliol-2021_NS 2.090031         NS   Juliol-2021
## Campus-Juliol-2021_S  1.806531          S   Juliol-2021
## Campus-Jun-2021_NS    2.701045         NS      Jun-2021
## Campus-Jun-2021_S     2.080060          S      Jun-2021
## Campus-maig-2020_NS   2.199104         NS     maig-2020
# Fer el model
model_shannonabellescampus <- lmer(Shannon ~ Tractament + (1|Data_Mostreig), data = taula_shannon_abellescampus)
summary(model_shannonabellescampus)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Tractament + (1 | Data_Mostreig)
##    Data: taula_shannon_abellescampus
## 
## REML criterion at convergence: 1.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.03570 -0.47827  0.03651  0.57598  0.80378 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  Data_Mostreig (Intercept) 0.03686  0.1920  
##  Residual                  0.02461  0.1569  
## Number of obs: 7, groups:  Data_Mostreig, 4
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   2.4004     0.1240  3.9717  19.363 4.42e-05 ***
## TractamentS  -0.5097     0.1248  2.2656  -4.083   0.0444 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.398

TractamentS: quan el tractament és S, s’espera que l’índex de Shannon disminueixi 0,5097 respecte el tractament NS. És un resultat significatiu! (p-valor= 0.0444).

# Plot del model
plot(fitted(model_shannonabellescampus), resid(model_shannonabellescampus), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

MODELS LINEALS D’ANÀLISIS FUNCIONALS

Inter-tegular span (ITS) (model lineal mixt)

Abelles totals, amb Zona com a factor aleatori

#Obrir la base de dades amb la columna de les mesures de la ITS i filtrar dades per abelles
abelles_ITS <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
abelles_ITS <- abelles_ITS %>%
  mutate(data_mostreig = paste(Mes, Any, sep = "-"))

head(abelles_ITS)
## # A tibble: 6 × 27
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 19 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>, data_mostreig <chr>
# Model lineal mixt ITS (mida)
abelles_ITS_filtrades <- abelles_ITS[!is.na(abelles_ITS$Inter_tegular_span),]
model_abelles_ITS <- lmer(`Inter_tegular_span` ~ Tractament + (1|Zona) + (1 | data_mostreig), data = abelles_ITS_filtrades)
summary(model_abelles_ITS)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | Zona) + (1 | data_mostreig)
##    Data: abelles_ITS_filtrades
## 
## REML criterion at convergence: 1425.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3887 -0.6319 -0.1823  0.1260  5.6755 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 0.03817  0.1954  
##  Zona          (Intercept) 0.01395  0.1181  
##  Residual                  0.45922  0.6777  
## Number of obs: 683, groups:  data_mostreig, 6; Zona, 3
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   1.69083    0.12062   1.57957  14.018    0.012 *  
## TractamentS  -0.40651    0.06206 592.97166  -6.551 1.24e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.230

La ITS de les abelles disminueix -0.40651 quan el tractament és S (abelles més petites quan està segat). És significatiu (p-valor = 1.24e-10)

# Plot del model
plot(fitted(model_abelles_ITS), resid(model_abelles_ITS), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Visualitzar gràficament amb boxplot

ggplot(abelles_ITS_filtrades, aes(Tractament, Inter_tegular_span, fill = Tractament)) +
  geom_boxplot() +
  facet_wrap(~Zona) +
  labs(x = "Tractament", y = "ITS") +
  scale_fill_manual(values = c("S" = "#DB6443", "NS" = "#74AA90")) +
  theme_minimal()

Abelles del Campus, amb data_mostreig com a factor aleatori

# Filtrar les dades només per abelles del Campus
abellescampus_ITS <- abelles_ITS[abelles_ITS$Zona == "Campus",]

#Definir nova variable: "data_mostreig", combinació de mes i any de mostreig
abellescampus_ITS <- abellescampus_ITS %>%
  mutate(data_mostreig = paste(Mes, Any, sep = "-"))

head(abellescampus_ITS)
## # A tibble: 6 × 27
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 19 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>, data_mostreig <chr>
# Fer el model
model_abellescampus_ITS <- lmer(`Inter_tegular_span` ~ Tractament + (1|data_mostreig), data = abellescampus_ITS)
summary(model_abellescampus_ITS)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | data_mostreig)
##    Data: abellescampus_ITS
## 
## REML criterion at convergence: 1196.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3101 -0.5834 -0.2132  0.1382  5.4645 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 0.04629  0.2151  
##  Residual                  0.42778  0.6540  
## Number of obs: 594, groups:  data_mostreig, 4
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   1.54523    0.11430   3.29968  13.519 0.000529 ***
## TractamentS  -0.36675    0.06534 537.59881  -5.613 3.19e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.243

Molt significatiu: la ITS de les abelles (mida) disminueix -0.36675 quan el tractament és S. Les abelles son més petites quan està segat.

# Plot del model
plot(fitted(model_abellescampus_ITS), resid(model_abellescampus_ITS), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

ITS abelles per espècies totals:

# Crear nova base de dades: només una espècie per cada tractament.
especies_ITS <- abelles_ITS %>%
  distinct(ID, Tractament, .keep_all = TRUE)
especies_ITS <- especies_ITS %>%
  mutate(data_mostreig = paste(Mes, Any, sep = "-"))
# Dialectus -- dialectus
head(especies_ITS)
## # A tibble: 6 × 27
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 19 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>, data_mostreig <chr>
# Fer el model (no podem considerar la zona com a factor aleatori ja que hem eliminat espècimens!)
model_especies_ITS <- lmer(`Inter_tegular_span` ~ Tractament +(1|Zona) + (1|data_mostreig), data = especies_ITS)
summary(model_especies_ITS)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | Zona) + (1 | data_mostreig)
##    Data: especies_ITS
## 
## REML criterion at convergence: 203.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3366 -0.6971 -0.2923  0.5604  3.1878 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 0.07046  0.2655  
##  Zona          (Intercept) 0.07477  0.2734  
##  Residual                  1.00207  1.0010  
## Number of obs: 70, groups:  data_mostreig, 6; Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)   2.2408     0.2593  2.1034   8.642   0.0112 *
## TractamentS  -0.1192     0.2619 65.7501  -0.455   0.6505  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.318

No significatiu

# Plot del model
plot(fitted(model_especies_ITS), resid(model_especies_ITS), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

ITS abelles per espècies del campus

# Crear nova base de dades: només una espècie per cada tractament.
especiescampus_ITS <- abellescampus_ITS %>%
  distinct(ID, Tractament, .keep_all = TRUE)
# Dialectus -- dialectus
head(especiescampus_ITS)
## # A tibble: 6 × 27
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 19 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>, data_mostreig <chr>
# Fer el model (no podem considerar la zona com a factor aleatori ja que hem eliminat espècimens!)
model_especiescampus_ITS <- lmer(`Inter_tegular_span` ~ Tractament + (1|data_mostreig), data = especiescampus_ITS)
summary(model_especiescampus_ITS)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | data_mostreig)
##    Data: especiescampus_ITS
## 
## REML criterion at convergence: 168.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4350 -0.7065 -0.2774  0.5608  2.6641 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 0.07751  0.2784  
##  Residual                  0.91037  0.9541  
## Number of obs: 60, groups:  data_mostreig, 4
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   2.0006     0.2138  5.6703   9.358 0.000117 ***
## TractamentS  -0.1076     0.2634 56.2365  -0.408 0.684543    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.465

No significatiu.

# Plot del model
plot(fitted(model_especiescampus_ITS), resid(model_especiescampus_ITS), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Substrat de nidificació (GLMM binomial)

Nest_type en funció del tractament de les abelles totals (zona i data_mostreig = factor aleatori)

abelles_ITS<-abelles_ITS[!is.na(abelles_ITS$Nesting_type),]

abelles_ITS["Nius"] = 0
abelles_ITS$Nius[abelles_ITS["Nesting_type"]=="S"]="S"
abelles_ITS$Nius[abelles_ITS["Nesting_type"]!="S"]="Not_S"

# Taula de contingència
taula_contingencia1 <- table(abelles_ITS$Tractament, abelles_ITS$Nius)
print(taula_contingencia1)
##     
##      Not_S   S
##   NS    71 324
##   S     14 291
# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia1)
ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(x = "Tractament", y = "Nº", fill = "Nesting_type") +
  ggtitle("Nest type totals") +
  scale_fill_manual(values = c("S" = "#3D405B", "Not_S" = "#F2CC8F")) +
  theme_minimal()

# Model:
abelles_ITS$Nius <- factor(abelles_ITS$Nius)
# Nomes Zona o tambe mostreig, com a factors aleatoris
model_abelles_nestt <- glmer(Nius ~ Tractament + (1|Zona) + (1|data_mostreig), data = abelles_ITS, family = binomial)
## boundary (singular) fit: see help('isSingular')
summary(model_abelles_nestt)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Nius ~ Tractament + (1 | Zona) + (1 | data_mostreig)
##    Data: abelles_ITS
## 
##      AIC      BIC   logLik deviance df.resid 
##    493.3    511.5   -242.6    485.3      696 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7820  0.2091  0.2282  0.4715  0.5146 
## 
## Random effects:
##  Groups        Name        Variance  Std.Dev. 
##  data_mostreig (Intercept) 3.795e-02 0.1947978
##  Zona          (Intercept) 1.277e-10 0.0000113
## Number of obs: 700, groups:  data_mostreig, 6; Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.4843     0.1645   9.024  < 2e-16 ***
## TractamentS   1.6262     0.3450   4.714 2.43e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.442
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Significatiu! Les abelles del tractament S tenen més probabilitat (1.6262) de fer nius al sòl que de fer nius a altres llocs que no siguin el sòl. És un resultat significatiu! (p-valor = 2.43e-06).

# Plot del model
plot(fitted(model_abelles_nestt), resid(model_abelles_nestt), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Nest_type en funció del tractament de les abelles del Campus (data_mostreig = factor aleatori)

abellescampus_ITS<-abellescampus_ITS[!is.na(abellescampus_ITS$Nesting_type),]

abellescampus_ITS["Nius"] = 0
abellescampus_ITS$Nius[abellescampus_ITS["Nesting_type"]=="S"]="S"
abellescampus_ITS$Nius[abellescampus_ITS["Nesting_type"]!="S"]="Not_S"

# Taula de contingència
taula_contingencia2 <- table(abellescampus_ITS$Tractament, abellescampus_ITS$Nius)
print(taula_contingencia2)
##     
##      Not_S   S
##   NS    63 291
##   S     13 243
# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia2)

ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(x = "Tractament", y = "Nº", fill = "Nesting_type") +
  ggtitle("Nest type Campus") +
   scale_fill_manual(values = c("S" = "#3D405B", "Not_S" = "#F2CC8F")) +
  theme_minimal()

# Model:
abellescampus_ITS$Nesting_type <- factor(abellescampus_ITS$Nius)
model_abellescampus_nest <- glmer(Nesting_type ~ Tractament + (1|data_mostreig), data = abellescampus_ITS, family = binomial)

summary(model_abellescampus_nest)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Nesting_type ~ Tractament + (1 | data_mostreig)
##    Data: abellescampus_ITS
## 
##      AIC      BIC   logLik deviance df.resid 
##    440.2    453.4   -217.1    434.2      607 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5162  0.2291  0.2360  0.4689  0.4997 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 0.02712  0.1647  
## Number of obs: 610, groups:  data_mostreig, 4
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.4968     0.1764   8.484  < 2e-16 ***
## TractamentS   1.5004     0.3860   3.887 0.000101 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.512

Significatiu! Les abelles del tractament S tenen més probabilitat (1.5004) de fer nius al sòl que de fer nius a altres llocs que no siguin el sòl. És un resultat significatiu! (p-valor = 0.000101).

# Plot del model
plot(fitted(model_abellescampus_nest), resid(model_abellescampus_nest), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Lèctia (GLMM binomial)

Lecty en funció del tractament de les abelles totals (Zona i data_mostreig = factors aleatoris)

# Taula de contingència
taula_contingencia3 <- table(abelles_ITS$Tractament, abelles_ITS$Lecty)
print(taula_contingencia3)
##     
##      Oligolectic Polylectic
##   NS         142        249
##   S           38        263
# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia3)

ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(x = "Tractament", y = "Nº", fill = "Lecty") +
  ggtitle("Lecty totals") +
   scale_fill_manual(values = c("Polylectic" = "#3D405B", "Oligolectic" = "#F2CC8F")) +
  theme_minimal()

# Model:
abelles_ITS$Lecty <- factor(abelles_ITS$Lecty)
model_abelles_lecty <- glmer(Lecty ~ Tractament + (1|Zona) + (1|data_mostreig), data = abelles_ITS, family = binomial)
## boundary (singular) fit: see help('isSingular')
summary(model_abelles_lecty)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Lecty ~ Tractament + (1 | Zona) + (1 | data_mostreig)
##    Data: abelles_ITS
## 
##      AIC      BIC   logLik deviance df.resid 
##    589.1    607.3   -290.6    581.1      688 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1853 -0.8175  0.1735  0.5793  1.2233 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 4.057    2.014   
##  Zona          (Intercept) 0.000    0.000   
## Number of obs: 692, groups:  data_mostreig, 6; Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   1.5679     0.8603   1.823   0.0684 .
## TractamentS   0.6398     0.2749   2.328   0.0199 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.119
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Les abelles tenen una major probabilitat de ser polilèctiques (0,6398) quan estan sotmeses al tractament S. És un resultat significatiu (p-valor= 0,0199), la qual cosa suggereix un efecte del tractament en la lecticitat de les abelles.

# Plot del model
plot(fitted(model_abelles_lecty), resid(model_abelles_lecty), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Lecty en funció del tractament de les abelles del campus (data_mostreig = factor aleatori)

# Taula de contingència
taula_contingencia4 <- table(abellescampus_ITS$Tractament, abellescampus_ITS$Lecty)
print(taula_contingencia4)
##     
##      Oligolectic Polylectic
##   NS         127        223
##   S           20        232
# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia4)

ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(x = "Tractament", y = "Nº", fill = "Lecty") +
  ggtitle("Lecty Campus") +
  scale_fill_manual(values = c("Polylectic" = "#3D405B", "Oligolectic" = "#F2CC8F")) +
  theme_minimal()

# Model:
abellescampus_ITS$Lecty <- factor(abellescampus_ITS$Lecty)
model_abellescampus_lecty <- glmer(Lecty ~ Tractament + (1|data_mostreig), data = abellescampus_ITS, family = binomial)

summary(model_abellescampus_lecty)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Lecty ~ Tractament + (1 | data_mostreig)
##    Data: abellescampus_ITS
## 
##      AIC      BIC   logLik deviance df.resid 
##    468.8    482.0   -231.4    462.8      599 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2618  0.0518  0.1706  0.5796  1.2245 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 6.291    2.508   
## Number of obs: 602, groups:  data_mostreig, 4
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   2.2170     1.3182   1.682   0.0926 .
## TractamentS   0.6375     0.3271   1.949   0.0513 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.091

Les abelles tenen una major probabilitat de ser polilèctiques quan estan sotmeses al tractament S. És un resultat marginalment significatiu (p-valor = 0.0513), la qual cosa suggereix un cert efecte del tractament en la lecticitat de les abelles.

# Plot del model
plot(fitted(model_abellescampus_lecty), resid(model_abellescampus_lecty), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")